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ABSTRACT 


A numerical  method  is  developed  to  compute  the  pressure  distribution  and 
normal  approach  in  a generalized  elliptical  contact  between  layered  linearly 
elastic  solids.  The  computed  quantities  are  obtained  as  an  approximate 
solution  to  an  integral  equation  formulated  for  the  most  general  case  of 
Hertz's  assumptions,  i.e.,  the  frictionless  contact  between  arbitrary  sur- 
faces whose  undeformed  normal  separation  can  be  approximated  by  the  separa- 
tion between  an  elliptical  paraboloid  and  the  tangent  plane  at  its  vertex. 

The  numerical  method  is  based  on  a discretized  representation  of  the  unknown 
pressure  distribution.  Values  of  pressure  are  defined  at  points  on  a two- 
dimensional  rectilinear  grid  and  a linear  function  approximates  the  pressure 
distribution  between  adjacent  points.  Expression  of  the  integral  equation 
at  each  grid  point  yields  a system  of  linear  equations  in  the  discrete 
pressures,  with  coefficients  formed  from  Integrals  of  the  kernel  function 
over  right  triangular  regions  between  adjacent  points.  The  kernel,  given 
by  a Hankel  transform  inversion,  is  integrated  numerically  by  Gaussian 
quadrature . 

The  method  is  applied  to  the  contact  between  a homogeneous  solid  and  a layered 
solid  consisting  of  an  isotropic  surface  layer  of  uniform  thickness  in  per- 
fect adheslan  to  an  isotropic  substrate.  Comparisons  with  available  solutions 
for  the  limiting  cases  of  Hertz  contact  between  homogeneous  solids  and  axi- 
symmetrlc  contact  of  layered  solids  establish  the  accuracy  of  the  numerical 
method.  Solutions  obtained  for  tliree  sets  of  material  properties  over 
ranges  of  the  contact  zone  ellipticlty  ratio  and  major  halfwidth  to  layer 

xi 


thickness  ratio  are  given  as  plots  which  show  the  najor  and  minor  halfwidths, 


approach,  and  central  pressure  as  functions  of  layer  thickness  at  fixed  load 
and  fixed  principle  curvature  ratio  of  the  parabolic  separation.  As  the 
layer  thickness  decreases,  the  contact  zone  ellipticity  ratio  is  found  to 
increase  for  a high-modulus  layer  bonded  to  a low-modulus  substrate  and 
decrease  for  a low-modulus  layer  bonded  to  a high-modulus  substrate. 
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a = major  contact  zone-  ha  If  width  (in.). 

a = tiwajor  contact  zone  halfwidth  for  axisvmmetric  Hertz  contact,  defined 
in  Kq nation  3.27  (in.). 

b = minor  contact  zone  lialfwidth  (in.). 

d = normal  approacii  (in.). 

d = normal  approach  for  axisymmeiric  Hertz  contact,  defined  in 
Equation  3.2  3 (in.). 

2 

Ej  = elastic  modulus  of  layer  (Ih/in.  ). 

2 

E^  = elastic  modulus  of  substrate  (Ib/in.  ). 

2 

E^  = elastic  modulus  of  homogeneous  solid  (Ib/ln.  ). 

F = magnitude  of  concentrated  force  (lb). 

g,  (i  = 1,  ....  6)  = constants  for  Hankel  transform  solution,  defined 

in  Table  2-1. 

G = Hankel  transform  of  stress  function,  defined  in  Equation  2.6. 
h = layer  thickness  (in.). 

.1  = Bessel  function  of  the  first  kind  of  order  zero, 

o 

= Bessel  function  of  the  first  kind  of  order  one. 

m(i)  (i  = 1,  ...,  n)  = grid  point  number  in  q direction  locating  contact 

zone  boundary. 

n = number  of  grid  points  along  major  halfwidth. 

N = total  number  of  grid  points  in  contact  zone. 

2 

p = contact  pressure  (Ib/ln.  ) . 

p^  = pressure  at  center  of  contact  (Ib/in.^). 

o pressure  at  center  of  contact  for  axlsymmetrlc  Hertz  contact,  defined 
^ o Eqtiation  3.26  (Ib/ln^). 

r = radial  coordinate  (in.). 

* major  radius  of  curvature  of  separation  profile  (in.). 


In 


xlli 


= minor  radius  of  ourvatiire  of  separation  prt)filf  (in.). 


i; 

y 

li^  = radial  (tangential)  displacement  (in.). 

11  = axial  (normal)  displacement  (in.), 

z 

u = normal  displacement  at  surface  (in.). 

" 2 
u = dimensionless  surface  displacement  (=  2R  u /a  ). 
o X o 

W = contact  load  (lb). 

W = dimensionless  contact  load,  defined  in  Equation 
X = coordinate  parallel  to  major  halfwidth  (in.). 

X = Love's  stress  function. 

X = modified  stress  function  (=  X/h^). 
y = coordinate  parallel  to  minor  halfwidth  (in.), 
z = axial  (normal)  coordinate 

a = major  halfwidth  to  layer  thickness  ratio  (-  a/h) . 

S = elastic  property  ratio  between  layer  and  homogeneous  solid,  defined 
in  Equation  3.6. 

r = dimensionless  normal  approach,  defined  in  Equation  3.9. 

£ = curvature  ratio  (=  R /R  )'. 

X y 

5 = dimensionless  coordinate  (=  z/h) . 
ri  = dimensionless  coordinate  (=  y/a). 

X = d imensi onless  coordinate  (=  r/a) . 

Vj  = Poisson's  ratio  of  layer. 

= Poisson's  ratio  of  substrate. 

= Poisson's  ratio  of  homogeneous  solid. 

^ = dimensionless  coordinate  (=  x/a). 
p =»  ellipticity  ratio  (=  a/b)  . 

p = dimensionless  coordinate  (=  r/h)  . T 

2 

T » radial  normal  stress  (Ib/in.  ). 
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= axial  normal 
= shear  stress 
= circumferent 
d imens  ionless 
dimensionless 


2 

stress  ( 1 b/ in . ) . 

(Ib/in.^) . 

2 

ial  normal  stress  (Ib/in.  ). 

pressure,  defined  in  Equation 

pressure  at  center  of  contact 

n;  j = 1,  m(i))  = values 

at  grid 


3.9. 


of  dimensionless 
points . 


pressure 


w = 

Hankel 

transform 

variable . 

= grid 

spacing 

in 

^ direction 

(=l/(n-l)) 

An 

= grid 

spac ing 

i n 

n direction 

(=Ar,/p)  . 
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PART  1 


INTR.ODUCrj^ 

Hertz's  classical  solution  (l-3|*for  the  contact  stresses  between  de- 
formable solids  is  frequently  applied  in  the  engineering  analysis  of  load 
bearing  mechanical  elements  such  as  ball  and  roller  bearings,  cams,  and 
gears.  Under  Hertz's  basic  assumptions,  the  equations  of  linear  elasticity 
apply,  friction  is  absent,  and  each  solid  is  treated  as  a halfspace 
under  normal  surface  loading.  In  the  contact  zone,  the  relative 
normal  surface  displacement  between  the  solids  is  prescribed  by  their 
deformation  into  contact.  The  contact  pressure  distribution,  which 
is  unknown,  depends  on  this  relative  displacement  through  a singular 
integral  equation  whose  kernel  is  the  surface  displacement  of  a half- 
space due  to  a concentrated  normal  force.  When  the  contact  is  between 
rounded  surfaces  over  a small  enough  region,  the  undeformed  normal 
separation  can  be  approximated  by  a quadratic  function;  i.e.,  the 
separation  is  that  between  a plane  and  an  elliptical  paraboloid.  The 
relative  displacement  profile  in  the  contact  zone  is  thus  parabolic. 

For  homogeneous  isotropic  solids,  the  pressure  profile  that  satisfies 
the  resulting  integral  equation  was  shown  by  Hertz  [1-3]  to  be  a 
semi-ellipsoid  whose  major  to  minor  axis  ratio  in  the  surface  of 
contact  depends  purely  on  the  shape  of  the  elliptical  paraboloid 
defining  the  undeformed  separation.  The  problem  may  be  formulated 
in  the  same  manner  when  one  of  the  solids  consists  of  an  isotropic 
surface  layer  backed  by  an  Isotropic  substrate  of  a different  material. 
However,  the  integral  equation  will  not  have  a simple  analytical  solu- 
tion because  the  concentrated  force  solution  for  a layered  halfspace 
is  given  by  a Fourier  or  Hankel  integral  that  requires  numerical  evalua- 
tion. 

Current  interest  in  the  contact  problem  for  layered  solids  follows 
from  a need  for  analytical  prediction  of  stresses  and  displacements  in 
coated  bearing  elements.  In  some  applications,  a hard  (high  modulus) 
coating  such  as  a ceramic  may  be  used  to  increase  the  surface  durability 
of  an  element;  in  others,  a soft  (low  modulus)  coating  such  as  an 

^Numbers  in  brackets  denote  references  listed  at  the  end  of  the  report. 
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elastomer  can  protect  elements  that  experience  impact.  Water-lubricated 
marine  bearings  and  gas-lubricated  directional  control  gyroscopes  are  sub- 
jects of  current  research  that  fall  into  the  latter  category.  Both  appli- 
cations involve  a bearing  surface  consisting  of  a thin  compliant  elastomeric 
layer  bonded  to  a relatively  rigid  substrate.  Contact  pressure  distributions 
must  be  determined  before  one  can  compute  the  subsurface  stresses  developed 
in  the  layer  under  concentrated  contact.  When  there  is  a large  enough 
difference  in  elastic  properties  between  layer  and  substrate,  the  true 
contact  pressures  can  differ  significantly  from  those  predicted  by  the 
Hertz  solution. 

To  date,  any  layered  contact  solutions  available  in  the  literature 
are  limited  to  plane  strain  and  axisymmetry  [4-8].  For  plane  strain, 
the  undeformed  separation  is  that  between  a plane  and  an  infinitely 
long  parabolic  cylinder;  the  kernel  is  the  solution  for  a line  load. 

For  axisymmetry,  the  undeformed  separation  is  that  between  a plane  and 
a paraboloid  of  revolution;  the  kernel  is  the  solution  for  a ring 
load.  Except  for  differences  in  the  mathematical  form  of  the  kernel, 
the  integral  equations  for  the  two  cases  are  the  same.  The  unknown 
pressure  is  one-dimensional  and  the  same  general  strategies  can  be 
used  to  obtain  approximate  numerical  solutions  to  both  problems.  There 
are  no  solutions  in  the  literature  for  the  general  case  of  three-dimensional 
contact  for  layered  solids,  in  which  the  undeformed  separation  is  that 
between  a plane  and  an  elliptical  paraboloid.  Here  the  kernel  is  the 
solution  for  a point  load,  and  the  unknown  pressure  is  two-dimensional. 

The  aim  of  the  present  study  is  to  solve  this  problem  numerically  by 
extending  to  two  dimensions  a method  that  has  been  used  for  plane  and 
axisymmetric  problems. 

The  literature  up  to  1972  on  plane  and  axisymmetric  layered  contact  is 
listed  by  Chen  and  Engel  [4].  Of  the  works  listed,  only  that  of  Tu  [5] 
has  direct  bearing  on  the  method  used  in  the  present  study.  In  obtaining 
a numerical  solution  to  the  axisymmetric  contact  of  a plate  pressed 
between  two  spheres,  Tu  represents  the  contact  pressure  distribution  by 
values  defined  on  a discrete  set  of  concentric  circles  with  the  pres- 
sure varying  linearly  in  radius  between  adjacent  circles.  The  analogous 
treatment  for  plane  problems  was  recently  put  forth  by  Gupta  and  Walowit 
[6]  who  treat  the  plane  strain  contact  between  a homogeneous  elastic 
solid  and  a layered  elastic  solid.  The  plane  contact  pressure  distribution 
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is  represented  by  values  defined  at  a discrete  set  of  points  with  t iu;  pres- 
sure varying  linearly  between  adjacent  points.  In  both  cases,  tiie  exact 
displacement  at  each  discrete  point  due  to  the  approximate  pressure 
distribution  is  obtained  in  terms  of  Fourier  integrals  and  matched  to 
the  prescribed  contact  displacement  to  produce  a set  of  linear  equations 
in  the  unknown  discrete  pressure  values. 

Chen  and  Engel  (A),  in  treating  the  axlsymmetric  contact  between  a 
rigid  paraboloid  and  a layered  elastic  solid,  represent  the  unknown 
pressure  distribution  p(r)  as  a truncated  series  of  base  pressures: 
p^(r)  = A^(l-r^/a^)^  (where  are  undetermined  coefficients  and 

a is  the  contact  radius).  With  exact  surface  displacements  due  to 

each  base  function  computed  from  Fourier  integrals,  the  coefficients 
A^  are  determined  by  matching  computed  displacements  to  prescribed 

contact  displacements  according  to  an  integral-least  square  or  collo- 
cation criterion.  Methods  similar  to  this  one  are  applied  in  many  of 
the  earlier  works  cited  by  Chen  and  Engel.  In  some  instances,  the 
form  of  the  base  pressures  follows  from  a power  series  approximation 
of  the  kernel  function.  Here  the  calculated  surface  displacement  due 
to  the  approximate  pressure  distribution  is  also  approximate.  Wu  and 
Ling  [7]  use  this  approach  to  obtain  approximate  solutions  for  the 
axlsymmetric  contact  between  a rigid  paraboloid  and  an  elastic  plate 
resting  on  a rigid  foundation. 

Methods  that  involve  series  approximations  of  the  kernel  are  difficult 
to  implement  when  the  layer  is  very  thin  relative  to  the  contact  zone 
dimension.  The  number  of  terms  needed  for  convergence  Increases  rapidly 
as  the  layer  thickness  decreases  [5].  Also,  in  the  case  of  a hard  layer 
on  a soft  substrate,  the  pressure  profile  shape  can  depart  significantly 
from  the  elllotical  shape  of  the  Hertz  solution.  Pressure  profiles  given 
in  References  [A]  and  [6]  for  this  case  show  a peak  near  the  edge  of 
contact  that  becomes  sharper  as  the  layer  thickness  decreases.  An 
unwieldy  number  of  base  functions  is  presumably  required  to  represent 
such  a profile  accurately.  The  discrete  pressure  methods  of  Tu  [5]  and 
of  Gupta  and  Walowlt  [6|  do  not  suffer  from  these  limitations  though  they 
do  require  more  computational  labor  than  metht)ds  based  on  scries  approxi- 
mations. The  kernel  evaluation  is  exact  and  pressure  Irregularities  can 
be  handled  by  increasing  the  number  of  discrete  points  in  troublesome  rc-glons. 
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In  the  present  study,  a natural  extension  of  the  above  discrete  pressure 
methods  is  applied  to  the  general  three-dimensional  contact  between  a 
layered  elastic  solid  and  a homogeneous  elastic  solid.  The  layered 
solid  consists  of  a single  isotropic  layer  in  perfect  adherence  to  an 
isotropic  substrate.  Since  the  unknown  pressure  distribution  depends 
on  two  space  variables,  it  is  represented  by  values  defined  at  discrete 
points  on  a two-dimensional  grid.  Following  Tu  [5]  and  Gupta  and 
Walowlt  [6],  adjacent  grid  points  delineate  elements  of  linearly  varying 
pressure,  but  this  time  in  two  dimensions.  The  basic  contact  zone 
division  is  a right  triangle;  the  elemental  pressure  distribution  is 

a plane.  The  point  load  solution  for  the  kernel  comes  from  a Hankel 

transform  analysis  similar  to  that  of  Lamb,  Terezawa,  and  Sneddon  [8]. 

As  in  the  one-dimensional  methods,  the  exact  evaluation  of  the  kernel 

makes  possible  the  exact  computation  of  surface  displacements  due  to 
the  approximate  pressure  distribution.  The  procedure  of  matching  com- 
puted displacements  to  prescribed  contact  displacements  at  each  grid 

point  is  then  employed  to  get  the  desired  system  of  linear  equations 

in  the  unknown  pressures.  The  contact  zone  boundary  prescribed  in 
setting  up  the  problem  is  approximate  since  by  definition  it  follows 
the  grid  lines  that  bound  the  outermost  pressure  elements.  Moreover, 
the  true  boundary  is  an  unknown  in  the  problem  because  it  depends  on 
the  mathematical  form  of  the  pressure  distribution.  Since  the  approxi- 
mate boundary  has  to  be  specified  to  calculate  the  approximate  pressure 
distribution,  the  best  approximate  boundary  for  a given  grid  size  is 
determined  by  an  iterative  process. 

Numerical  solutions  are  presented  for  three  distinct  material  combinations 
over  ranges  of  the  contact  zone  halfwidth  to  layer  thickness  ratio  and 
the  contact  zone  aspect  or  ellipticity  ratio.  A discrete  boundary  pre- 
scribed to  represent  an  ellipse  gave  satisfactory  results  in  all  cases. 
Checks  on  the  accuracy  of  the  method  include  comparison  of  computed 
results  with  the  Hertz  solution  [1-3]  and  numerical  results  of  Chen  and 
Engel  (4).  The  solutions  are  presented  graphically  In  dimensionless 
form  with  major  and  minor  halfwidths,  relative  appr:»ach,  and  central 
pressure  shown  at  a fixed  load  as  functions  of  layer  thickness  and 
ratio  of  the  principal  radii  of  curvature  of  the  parabolic  separation 
prof ile . 
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The  analytical  foundations  of  the  study  are  presented  in  Parts  2 and  i. 
Part  2 covers  the  Hankel  transform  analysis  leading  to  the  point  load 
solution;  Part  3 gives  the  derivation  of  the  integral  equation  in  dimen- 
sionless form.  Details  of  the  numerical  solution  procedure  for  the 
integral  equation  are  given  in  Part  4,  and  Part  5 contains  the  numerical 
results.  A simple  example  at  the  end  of  Part  5 demonstrates  the  appli- 
cation of  the  numerical  results  to  an  engineering  problem. 
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PARX_2 

SOLUTION  FOR  SURFACE  DISPLACEMENT  DUE  TO  CONCENTRATED  NORMAL  FORCE 


2.1  Formulation  of  Boundary  Value  Problem 

Consider  the  layered  elastic  halfspace  shown  in  Figure  2-1  consisting  of  an 
isotropic  layer  of  thickness  h with  Young's  Modulus  and  Poisson's  Ratio 
bonded  uniformly  to  an  infinite  isotropic  substrate  of  properties  The 

halfspace  is  described  by  the  Cartesian  coordinates  x,  y,  z,  and  cylindrical 
coordinates  r,  6,  z,  with  the  x-y  and  r-6  planes  coinciding  with  the  free  sur- 
face, and  with  the  z axis  directed  into  the  halfspace.  In  this  section,  an 
expression  will  be  derived  for  the  normal  displacement  at  any  point  on  the 
surface  due  to  a concentrated  force  F acting  normally  inward  at  the  origin. 
This  solution  will  form  the  kernel  of  the  integral  equation  for  the  contact 
problem  to  be  dealt  with  subsequently. 


In  cylindrical  coordinates  r,  z,  begin  with  Love's  stress  function  X for 
axlsymmetric  problems  [2 ] satisfying 


4 2 2 

V X = V'^(V^X) 


= 0 
1 _3 
r 3r 


+ 


3z 


(2.1) 


with  stress  components  expressible  in  terms  of  X: 

2 


(vV‘ 


) X 


3r 


rr  3z 

3 , „2  13,^ 

■ 3z  ' r 3r  ^ ^ 


T = ^ 1(2-  v)V^  - -^  ] X 

3z 

^ [(1-  v)v2  - ^2  1 " 


"re  = "ez  = ° 


and  displacement  components: 
1 + ^ 


u = - 
r 


u = 
z 


3^X 


E 3r3z 
1 + V 


3^. 


[2(l-v)v‘  - ^ 1 X 


3z 


(2.2) 


(2.3) 
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It  can  be  verified  by  direct  substitution  that  Equations  2.3  satisfy  the  dis- 
placement form  of  the  static  equilibrium  equations  of  isotropic  linear  elasticity 
when  X satisfies  Equations  2.1.  By  successive  application  of  the  strain-displace- 
ment and  stress-strain  relations,  Equations  2.2  are  obtained  from  Equations  2.3. 


The  solution  for  a concentrated  normal  force  will  be  obtained  as  the  limit 
of  the  solution  for  an  arbitrary  axisymmetric  compressive  normal  loading  as 
the  area  over  which  it  acts  approaches  zero  while  the  total  load  remains  finite. 

For  the  analysis  which  follows,  it  will  be  convenient  to  work  with  the  coordinates 
z^,  originating  at  z = 0 and  z^  originating  at  z = h,  both  in  the  z direction. 

Also,  let  the  superscripts  (1)  and  (2)  denote  quantities  in  the  layer  and  substrate, 
respectively.  The  boundary  conditions  on  the  top  surface  for  a normal  loading 
p(r)  with  zero  shear  loading  can  now  be  expressed  as: 


"zz  ' ("l  = 


(1) 

:z 

(1) 


T ' ■ (z  = 0)  = 0 
rz  1 


For  the  case  of  perfect  bonding  which  this  analysis  will  be  concerned  with, 
the  continuity  conditions  at  the  layer-substrate  interface  are: 
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= 0) 
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^^1 
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= u(2) 

z 
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= 0) 

In  addition,  the  stress  components  must  satisfy 
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zz 

-> 
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as 

z + 00 
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rz 
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z,r  *■  «> 
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J*  ■>.  00 

(2.5) 


in  order  that  all  surface  tractions  vanish  at  “. 


2.2  Formal  Solution  by  Integral  Transforms 


In  many  cases,  partial  differential  equations  can  be  reduced  to  ordinary 
ones  through  the  use  of  integral  transforms.  For  axisymmetric  problems 
where  the  radial  variable  ranges  from  0 to  the  Hankel  Transform 
is  generally  applied.  The  formulation  used  here  is  similar  to  those 
of  Lamb,  Terezawa,  and  Sneddon  [ 8 ] in  obtaining  solutions  for  half- 
space and  plate  problems  with  axial  symmetry. 

Lt't  f(r)  be  a function  whicli  satisfies  the  sufficient  conditions  for 
representation  by  the  Hankel  Integral  Formula  [ 8 ].  Its  Hankel  Trans- 
form and  inverse  transform  are  given  by 

00 

F(w)  = /rf(r)J  ((jor)dr 
o 

0 

00 

f(r)  = /u)F(u>)J  (wr)dui 
o 

0 


Before  taking  the  Hankel  Transforms  of  Equations  2.1  through  2.5,  it 

will  be  helpful  to  rewrite  them  in  terms  of  the  dimensionless  coor- 

— '3 

dinates p = r/h,  z/h,  and  the  stress  function  X = X/h  which  has  the 

dimension  of  stress.  Equations  2.1  become 


h " 2 2" 

V X = y (V  X)  = 0 


v2  = -Jl.  + i A + _9_ 
— 2 - 2 

3p  p 9p  9^ 


with  the  stresses  and  displacement  from  Equations  2.2  and  2.3 

,2 

d , t 
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= —■  (vV^  - X 

""  9^  9^2 


3 / n2  12" 
^eo  - p Tp-^"" 


r.  ■ i -,;r 


(2.1a) 


(2.2a) 
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h(l+v) 

u ^ - 

r E 

^ ^ h(l+v) 

z E 

The  boundary  and  continuity  conditions  from  Equations  2.4  and  2.5  can 
be  written  as 

(2.4a) 

(2.5a) 

with  and  ^2  defined  analogously  to  z^  and  Z2. 

Denoting  the  Hankel  Transform  of  X(p,C)  by 
00 

G(w,C)  = /pX(p,5)J  (pu))dp 
o 
0 

and  taking  the  transform  of  Equations  2.1a,  one  obtains 


which  has  the  general  solution 

G = (A  + B4)e““^  + (C  + D4)e“^  (2.6) 

where  the  constant  coefficients.  A,  B,  C,  and  D,  will  generally  depend 
on  u).  Taking  the  Hankel  transforms  of  Equations  2.2a  and  2.3a,  and 
then  inverting,  one  obtains  the  following  expressions  for  the  stresses 
and  displacements  in  terms  of  G. 
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oo  oo 

dG  ^ , 2 dG  , . 1 , 2,  ,dG  , 

T = /U'J  (po))  [v — - + (l-v)ij)  — J du)  - = ; oj  J,(poj)— - doj 
rr  o , i , p 1 clc 

0 d c dr,  ^ o 


, ^ , d G 2 dG  , . ^ 1 , 2,  . dG  , 

T.  _ wJ  (pw)  [ r -U)  ldu  + =/u)J(pu))—  dio 

90  o ,,3  d?  p 1 d^ 

o d ? 0 


T = /wJ  (pco)  [(1-v) — ^ - (2-v)w  -T-  ] du 
zz  o , d d 5 

o d ^ 


(2.7) 


T = /u^J  (pw)  [v^— ^ + (l-v)a)^G  ] doj 
o ^ dc 


hd+v)  , 2,  , dG  , 

u = / u)  J (pio)  dw 

r E i d ^ 
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u = (p^)  [(i-2v)  ^ - 2(l-v)w^G  ] di 

z E o , ,2 

o d 5 


The  function  G will  take  the  following  forms  in  the  layer  and  substrate, 
respectively : 


(1)  “^2^ 

g"-  ^ = (A^  + + (Cj^  + D^C^e  . 0 £ 1 1 

(2)  ~‘^^2  “^2 
G^  ' = (A^  + B,^r,^)e  + (G^  + 02^2)6  , 0 1 ^2  - " 


(2.8) 


V ) 


In  order  for  the  tractions  to  vanish  as  " *^2  ~ *^2  ~ remaining 

six  constants  are  obtained  from  the  two  boundary  and  four  continuity 
conditions  in  Equations  2.4a  and  2.5a.  Using  Equations  2.7  to  express 
these  conditions  in  terms  of  G gives: 

[(1-Vi)  ^ - (2-v^)  1 0 = 

d;^  11 


(2.9) 


[V  ^ ^ + (1-v  )w  Gd>  ] = 0 
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1+v 


1+v 


IT 


dC 


jf  1(1-2''!)^^  - 2(l-vpA<l' 
2(1-»,)A<^> 


where  P(a))  is  the  Hankel  transform  of  the  normal  loading  p(p) 


P(a))  = /pp(p)  J (pw)  dp 
o 

Substituting  from  Equations  2.8  into  Equations  2.9  yields  six  simul- 
taneous equations  in  the  six  unknown  constants  , A2, 

and  B2.  It  is  convenient  to  express  the  unknown  constants  as  the 
quantities : 


gj^  = - A^  [oj  /P(o))] 

82  = - [(jO^/P(a))  ] 

g^  = - [u?/P(ii))] 

g^  = - Dj  [(o^/P((o)] 

8^  = - A2  |(.)^/P((.0| 

g^  = - B2  (oj^/PCu))  ] 


(2.10) 


(2.11) 
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which  turn  out  to  depend  only  on  u),  2-1 

gives  the  resulting  simultaneous  equations  in  g^^  through  g^  in  matrix 
form. 

At  this  point.  Equations  2.7,  2.8,  2.10,  and  2.11,  and  Table  2-1  could 
be  applied  to  compute  the  stresses  and  displacements  at  any  point  in 
the  halfspace  due  to  the  normal  surface  loading  p(p).  However,  the 
only  quantity  needed  to  formulate  the  contact  problem  is  the  normal 
surface  displacement  = 0)  which  will  be  denoted  by  u^. 

Substituting  from  the  first  of  Equations  2.8  into  the  sixth  of 
Equations  2.7,  setting  = 0,  and  using  Equations  2.11  and  the  second 
equation  in  Table  2-1  to  simplify  the  result,  the  following  expression 
is  obtained: 

2h(l-v  “ 

u (p)  = f [g.((D)-g,  (a))]P(u)J  (pw)doj  (2.12) 

O Z H O 

1 0 

To  obtain  the  specific  form  of  P(u))  for  a concentrated  normal 
force  F at  the  origin,  begin  with  the  following  loading: 

P _<  a/h 

p > a/h 

J (au/h) 

for  the  Hankel  transform  of  this  normal  loading.  The  limit  of  this 
expression  as  a ->■  0 is 

P((i))  = — ^ 

2nh 

the  Hankel  transform  of  a normal  loading  consisting  of  a concentrated 
force  of  magnitude  F at  the  origin.  Setting  F = 1,  write  Equation 
2.12  for  this  particular  P(w) . The  desired  end  result,  the  surface 
displacement  due  to  a unit  concentrated  force,  is  now  obtained. 

Denoting  it  by  6^,  one  has: 


F/ira 


p(p)  - 


Equation  2.10  gives 
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TABI.K  2-1 


SIMULTANEOUS  EQUATIONS  IN  CONSTANTS  TOR^HANKEL  TRANSFORM  SOI,UTr_ON 
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1-v  <» 

^ n E ~li  ^ [g2(w)-g^(u)  ] J^(ps)da) 


(2.13) 


To  obtain  algebraic  solutions  for  and  g^  from  the  6x6  system  in 
Table  2-1,  first  solve  the  third  and  fourth  equations  for  e'^g^  and 
e*^g,  and  use  the  resulting  expressions  to  eliminate  g^  and  g from 

O JO 

the  fifth  and  sixth  equations.  Then  solve  the  first  and  second 
equations  for  g^^  and  g^,  and  use  the  resulting  expressions  to  eliminate 
gj^  and  g^  from  the  modified  fifth  and  sixth  equations.  The  resultant 
form  of  the  fifth  and  sixth  equations  may  now  be  solved  for  g2  and 
g^  to  yield: 


g2'a))  - g^Ccj)  = 


^l(‘^21  ^22^  ^12^ 

''ll‘^22  ■ ^12*^21 


(2.1A) 


where  = 3 - 4v^  ~ Y + (1+Y)e^*^ 


'*'12  (1+Y)2we 


2a) 


= (l-Y)(3-4v^-2aj)  + 2y(2v2  - 2v^)  + (1+Y)e 


- 2 
^22  = 1 Y + e 
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k2  = (1-y)  + (1+Y)e^“ 
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[2  + 2Y(2v2n^+2v^n2+Av^-l)  - (l+y)  (Av j^-l-2aj) 


with 


ni  = 1 - 2^ 
' '"2 


"2  = 1- 


Y = yO-Av^) 


(2.16) 


It  is  easily  seen  that  for  large  w: 

g2(w)  - g^(“)  = 1 + 0(w^e  ^“) 


(2.17) 


This  result  permits  recasting  of  the  improper  integral  in  Equation 

2.13  to  facilitate  numerical  evaluation.  Add  and  subtract  J (pw) 

o 

inside  the  integral  and  recall  tlie  basic  result 


f J (pw)  du)  = =r 
o p 


.15) 
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to  obtain 


tie  h ^ [82^'-^^ 
1 L o 


g,  (oj)  - 1]  J (p(*j)  doj  + 
k o 


fl 


The  l/p"  term  is  the  surface  displacement  for  a homogeneous  half- 
space of  properties  E^,  v^,  and  the  integral  contains  the  effect  of 
layering.  (Note  from  Equations  2.15  through  2.17  that  - 

g^(w)  -1  = 0 for  all  CO  when  = 1 and  It  is  clear 

from  Equation  2.17  that  the  improper  integral  in  Equation  2.18 
converges  like  (o^e  ^ . Thus,  numerical  evaluation  to  high  accuracy 


is  possible  using  an  upper  limit  co  of  moderate  magnitude. 

o 

error  is  given  by: 

~ -2(0 

/ [g-Cu)  - g.  ((o)  - 1]  J (p(o)  d(o  = 0 (co  e °) 
(0/4  o o 
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The 


(2.18) 


(2.19) 


In  this  study,  (o  was  assigned  a value  of  10,  leaving  an  error  of  the 
-7  ° 

order  of  10 


In  terms  of  the  dimensionless  Cartesian  coordinates; 

C = x/a 
r,  = y/a 

J 

where  a will  be  a characteristic  length  for  the  contact  problem. 

Equation  2.14  may  be  expressed  as 

1-v  2 

« (C.n)  - aK'(uC,cin) 

O ^ 2^  3 

where 

K'(c,n)  » /(g2(4))  - g^((o)  - l]J^(,o/<^^  + n )d(0  + > - j-  - - -2  (2.22) 
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For  notational  simplicity,  define: 


= riK'(ar,r,n) 


(2.24) 


to  give  4 the 
o 


4 
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f orm : 
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K’(7,n) 


(2.2'j) 


Equation  2.25  gives  the  normal  surface  displacement  at  ^,n  due  to 
a unit  concentrated  normal  force  at  the  origin.  For  the  same  force 
positioned  at  C'»n’,  the  displacement  at  C.n  is  simply: 

1-^1^  - 

(2.26) 

o TiE^a 

This  relation  will  be  used  in  Section  3.1  to  obtain  the  displacement 
due  to  a non-axlsymmetr ic  pressure  distribution,  with  the  function  K' 
forming  a portion  of  the  kernal  of  the  integral  equation  for  the 
contact  problem. 


It  is  important  to  note  from  Equations  2.20  and  2.23  that  the  function 
K'  is  singular  when  its  arguments,  the  C and  q positions  of  the  point 
in  question  with  respect  to  the  position  of  the  concentrated  force, 
vanish.  The  singularity  is  isolated  to  the  Inverse  square  root  term 
which  was  subtracted  from  the  Hankel  inversion  integral  in  order  to 
form  an  integrand  which  would  die  out  exponentially  for  large  uj. 

The  integral,  as  a result,  is  fully  convergent  for  all  values  of 
the  arguments  of  K'.  Letting  R be  the  distance  between  the  con- 
centrated force  and  the  point  in  question,  the  singular  term  is  1/R. 
This  is  an  Integrable  singularity  in  the  sense  that  the  singularity 
vanishes  when  Equation  2.24  is  Integrated  to  obtain  the  displacement 
due  to  a continuous  pressure  distribution  over  a finite  region. 


f' 
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PART_3 

ANALYTICAL  FORMULATION  OF  I NTEGRAL  EQUATION  FOR  CONXACT_PJ^OBLEM 

3 . 1 Surface  Displacement  Due  to  Non-Ax Isymmetrlc  Pressure  Distribution 

To  formulate  the  contact  problem,  an  expression  is  needed  for  the 
normal  surface  displacements  of  a layered  halfspace  due  to  a pressure 
distribution  such  as  that  shown  in  Figure  3-1.  The  pressure  is 
symmetric  about  the  x and  y axes  and  is  bounded  by  a closed  curve  of 
the  same  symmetry  with  vertices  at  (+  a,0)  and  (o,+b).  The  pressure 
vanishes  everywhere  on  and  outside  the  bounding  curve. 

Using  the  dimensionless  coordinates  C = x/a  and  n = y/a  defined 
previously  and  letting  n = +f(C),  (with  f(-0  = f(C),  f(0)  = b/a, 
and  f (±1)  = 0),  denote  the  upper  and  lower  portions  of  the  bounding 
curve.  Equation  2.26  can  be  used  to  set-up  the  following  integral 
expression  for  the  normal  surface  displacement  u^(C,n)  due  to  the 
pressure  distribution  in  Figure  3-1: 

a(l-v  2)  1 f'r') 

(C.n)  = — ^ p(C',n’)K’(c-C',n-n’)dn'd^'  (3.1) 

’'^l  -1  -f(c') 

which  is  the  continuous  sum  over  the  pressure  distribution  of  the 

displacements  due  to  concentrated  forces  at  (C',n')  of  magnitude 
2 

p(^'«n')a  dn'dC'-  As  discussed  at  the  end  of  Section  2.2,  u^(C.h) 
is  bounded  and  continuous  everywhere  despite  the  singularity  of  K' 
at  (^,n)  = (^'.n').  Of  course  the  pressure  p(C'.n')  must  be  continuous 
and  sufficiently  smooth, 

3.2  Derivation  of  Integral  Equation 

Under  the  classical  Hertzian  assumptions  [1-31,  the  problem  of  contact 
between  two  elastic  solids  can  be  reduced  to  a consideration  of  the 
contact  between  an  elliptical  paraboloid  and  a hal f space . Figure  3-2 
shows  the  situation  of  interest  here  in  which  the  paraboloid  is  a 
homogeneous  solid  of  properties  E^,v,,  and  the  halfspace  is  the  layered 
solid  analyzed  in  Part  2.  The  load  direction,  parabolic  axis  of 
the  paraboloid,  and  z axis  of  the  layered  half space  are  coincident. 

The  X and  y axes  on  the  surface  of  the  layered  halfspace  coincide. 
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respectively,  with  the  major  and  minor  elliptical  axes  of  the 
paraboloid.  In  the  unstressed  state,  the  two  solids  are  in  point 
contact  at  the  origin  with  the  initial  normal  separation  given  by 


where  R and  R are  the  radii  of  curvature  at  the  vertex  of  the  un- 
X y 

deformed  paraboloid  in  the  x-z  and  y-z  planes,  respectively.  When 
the  solids  are  pressed  together  under  a load  W,  the  separation  at 
any  point  x,y  is  decreased  by  the  normal  approach  between  the  two 
solids  and  increased  by  the  sum  of  their  local  normally  inward  sur- 
face displacements.  Thus,  the  final  separation  under  load  is 


U,.  = U,  - d + u^  + u 
f i o o 


where  d is  the  normal  approach;  u and  u are  the  surface  dis- 

o o 

placements  of  the  half space  and  paraboloid,  respectively.  The  two 
solids  will  be  in  contact  over  a finite  region  centered  at  the  origin 
in  which  = 0,  forcing  the  surface  displacement  to  satisfy: 


I 

u 

o 


+ 


(3.2) 


The  load  W is  transmitted  through  the  contact  zone  by  a pressure 
distribution  p(x,y)  as  shown  in  Figure  3-1  which  acts  normally 
inward  on  the  surface  of  each  solid,  vanishing  on  and  outside  the 
contact  zone  boundary.  As  in  the  statement  of  the  classical  Hertz 
problem,  no  shear  stress  is  transmitted  between  the  contacting  solids. 


Equation  3.1  holds  exactly  as  written  for  where  p(^,n)  is  the 

contact  pressure  distribution;  n = +f(^)  is  the  contact  zone  boundary, 
and  a and  b are  the  major  and  .minor  contact  zone  halfwidths.  The 
corresponding  surface  displacement  of  the  homogeneous  solid  is  given 
by: 


il 

^‘o  = -.T-- 


/ 


/ 


-1 


P ( h ' , n ' ) 


(n-n') 


dn'  df/ 

f \ ^ 


(3.3) 
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Adding  Equation  3.3  to  Equation  3.1  for 


I II  '^l  ^ ^ 

u (C,n)  + u (f.,n)  = ^ f 

o o ttE^  ^ 


u^  yields: 

f(C') 

/ P(C' ,n')K(C-C' ,n-n') 
-f(C’) 


dn'dr/ 

(3.4) 


where  the  function  K is  defined  in  terms  of  a basic  kernal  function 

K(C,n): 


K(5/i  ) = 


/[g2(w)  - g^(w)  - l]J^(w 
o 


n )dw  + 


7^ 


1 + 


with 


according  to: 


(3.5) 


(3.6) 


K(C,h)  = aK(aC,an)  (3.7) 

Approximating  the  infinite  integral  in  w as  discussed  in  Section  2.2, 
the  form  of  K(C,h)  to  be  used  in  the  subsequent  numerical  analysis  is: 

U)  _____ 

K(C.ri)  [g  (u))  - g (oi)  - 1]J  (a)/c^  + q^)da)  + 7LI-L  (3.8) 

o /2  ^ 2 

✓ S + h 

2 2co 

with  an  error  of  order  e ° associated  with  this  approximation. 

By  combining  Equation  3.1  with  Equation  3.4  and  defining  the  dimension- 
less quantities: 

2 


r = 2R  d/a 

X 

e = R /R 
X y 


(3.9) 


<)>  = 


2R^(l-v^")p 

TTE^a 
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one  obtains  the  integral  equation 
o 0 1 f(C') 

r-C-en  = / f <))(C',n’)K(c-C',n-n’)dn'dC'  (3.10) 

-1  -f(C’) 

which  holds  for  all  ^ and  n within  the  contact  zone.  Outside  the 
contact  zone  and  along  its  boundary,  ifCC.n)  = 0. 


Taking  advantage  of  the  symmetry  of  the  problem,  Equation  3.10  can 
be  written  as  an  integral  over  the  first  quadrant  of  the  contact 
zone : 

2 2 ^ 

r - C - en  = / / <^(C',n')Q(C,n;C',n')dn’dC'  (3.11) 

o o 

where 


Q(C,n;C',n')  = K (^-^',n-n')  + k (C  + C',n-n') 
+ K (C-^' .1-1  + n')  +K  (C+C',n+n') 


(3.12) 


A dimensionless  contact  load  W may  be  obtained  from  the  Integral  of 
())(C,ri)  over  the  contact  zone: 

1 f(C) 

W = A / / ())(C^, ri)dndC  (3.13) 

o o 

Taking  the  expression  for  the  actual  contact  load  W: 
a af(x/a) 

W = 4 / / p(x,y)dydx 

o o 

and  writing  p,  x,  and  y in  terms  of  (j>>  C.  and  n,  respectively,  gives 
the  definition  of  W as: 

_ 2R  (1-v  ^)W 

W = i (3.14) 

TTE^a"^ 

This  study  will  be  concerned  with  numerically  solving  a discretized 
form  of  Equation  3.11  to  obtain  approximate  distributions  of  ({i(C.h)  along 
with  r,  e,  and  W for  various  layered  contact  configurations.  Before 


w 


proceeding  with  the  numerical  formulation,  it  is  of  some  value  to 
consider  the  classical  Hertzian  solution  as  a closed  form  analytical 
solution  to  Equation  3.11  and  making  use  of  that  solution  to  define 
some  dimensionless  variables  more  amenable  to  direct  physical 
interpretation  than  cfi,  W,  and  F. 


The  problem  at  hand  reduces  to  the  classical  Hertzian  contact  problem 
for  homogeneous  elastic  solids  when  one  sets  E^/E^  = 1 and 
in  Equations  2.1A  through  2.16  for  [g2(w)  - g^((D)].  It  is  easily 
verified  that  under  these  conditions  g2(w)  - &^(w)  = 1 for  all  u)  , so 
that  the  Hankel  inversion  integral  in  Equation  3.5  for  K(?,q)  vanishes 
identically,  causing  Equation  3.10  to  reduce  to: 


r 


2 

en 


J (l+B)  ()>(€' .n')dn’dc’ 

+ (n-n')^ 


(3.15) 


This  is  the  integral  equation  solved  analytically  by  Hertz  [ 1 ] , 
with  the  well  known  result  consisting  of  an  ellipsoidal  pressure  dis- 
tribution and  elliptical  contact  zone  boundary: 


/2  2 2 
1 - C - P n 

f(0  = ^ v/l  - 


(3.16) 


where 


p i a/b 


(3.17) 


Due  to  the  ellipsoidal  form  of 
, o 


W = ^ 
3p 


(3.18) 


The  constants  e,  F,  and  41°  are  given  in  terms  of  complete  elliptic 
integrals  whose  arguments  depend  on  p: 


P^E(p)  - K(m) 
K(h)  - E(m) 


(3.19) 


^ I 


i 


I 


2 j 


t 


I 


I 


I' 

( ■ 

I 

1' 


, ^ MK(y) 

K(y)  - E(m) 

o ^ py 

Ti(H-B)[K(y)  - E(y)] 


where 


p>  1 


K(y)  and  E(y)  are  complete  elliptic  integrals  of  the  first  and 
second  kind,  respectively,  defined  by: 


2 -1/2 

K(y)  = / (1-ysin  9)  ' dG 


2 1/2 

E(y)  = ; (1-ysin  0)  ' d6 


(3.20) 


(3.21) 


(3.22) 


(3.23) 


Tables  and  polynomial  approximations  of  these  functions  are  found 
in  Reference  [ 9]. 

For  the  axisymmetric  case,  the  contact  zone  is  circular  and  p = 1. 

The  constants  reduce  to: 

e = 1 

r = 2 

.o  4 (3.24) 

<!>  = ^ 

1T^(l+3) 

W = — ■ 

3Tr(l+B) 

In  most  physical  situations  to  which  contact  theory  can  usefully  be 
applied,  the  load  W and  curvature  ratio  e are  known.  Typically,  one 
wishes  to  determine  the  normal  pressures  p(x,y),  major  and  minor  half- 
widths a,b,  and  approach  d.  Denoting  by  p^(x,y),  a^,  and  d^. 


the 
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values  of  these  quantities  for  the  case  of  an  axisymmetr ic 
Hertzian  contact,  where  the  layered  halfspace  shown  in  Figure  3-2  is 
homogeneous  in  the  layer  properties  the  second,  third,  and 

fourth  of  Equation  3.2A  may  be  equated,  respectively,  to  the  first 
and  third  of  Equation  3.9  and  Equation  3.1A  to  give; 


(3.25) 


Po  i p^(0,0)  = 


2a  E, 
o 1 


TlR^(l+g)(l-V^  ) 


■3R^W(H-e)(l-v^^)1 

^o  = L J 


1/3 


(3.26) 


(3.27) 


These  three  equations  can  be  combined,  respectively,  with  the  first 
and  third  of  Equation  3.9  and  Equation  3.1A  to  yield  the  dimension- 
less quantities: 


1 2 

d/d  = -X  (a/a  ) r 

0^0 

(3.28) 

p/Pq  " 1-  7r^(l+6)  (a/a^)<j> 

(3.29) 

a/a^  = 2[3it(H-6)w] 

(3.30) 

In  addition,  a dimensionless  minor  halfwidth  and  layer  thickness 
are  clearly  given  by: 

a/a 

b/a  = - (3.31) 

o p 


a /h 
o 


(3.32) 


Through  the  direct  mathematical  solution  to  Equations  3.11  and  3.13,  one 

obtains  the  dependent  variables  of  the  problem:  4i(C,ri,  W,  T,  and  e, 

as  functions  of  the  independent  variables:  Ej^/E^,  V2. 

V-,  a,  and  p.  However,  since  the  quantities  a , d , and  p°  are  functions 
3 ^ o o *^0 

of  known  physical  variables  through  Equations  3.25  through  3.27,  Equations 
3.28  through  3.32  provide  a useful  physical  interpretation  of  this 
solution,  whereby  the  dependent  variables  are  taken  as  p/p°« 


27 


\ 

I 


b/a  , and  d/d  , with  the  independent  variables  E,/E„,  E,/E,,  v , v 
o o 12  13  1 

v_,  a /h,  and  c. 

3 o 


When  the  solution  is  viewed  in  this  manner,  all  of  the  independent 
variables  are  functions  of  generally  known  physical  parameters: 

W,  h,  R^,  R^,  Ej^,  E2,  E^,  and  v^.  The  dependent  variables 

consist  of  the  unknowns:  p,  a,  b,  and  d,  each  scaled  on  a function 
of  the  known  physical  parameters.  Unfortunately,  a^/h  arid  e, 
while  Independent  variables  physically,  are  dependent  mathematically. 
Thus,  to  obtain  a solution  for  particular  values  of  a^/h  and  e, 
one  must  perform  an  interpolation  of  solutions  to  Equation  3.11  over 
ranges  of  the  mathematically  Independent  variables  a and  p. 


Equations  3.18  through  3.21  may  be  combined  with  Equations  3.28 
through  3.31  to  yield  for  elliptical  Hertzian  contact: 


(3.33) 


For  this  case.  Equation  3.19,  together  with  Equations  3.22  and  3.33, 
must  be  solved  iteratively  to  determine  p given  a value  of  e.  The 
desired  quantities  are  then  determined  directly  from  Equations  3.33 
and  3.31. 


I 


f 


r 
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?AiT  A 

NUMERICAL  FORMULATION 

4.1  Discretized  Integral  Equation  and  Method  of  Solution 


Equation  3.11,  which  is  the  integral  equation  to  be  solved  for  the 
contact  pressure  distribution  4>(^,n),  is  discretized  by  dividing 
the  first  quadrant  of  the  contact  zone  into  a rectangular  mesh  as 
shown  in  Figure  4-1  with  n-1  uniform  divisions  between  i = 0 and 
C = 1,  and  the  same  number  of  uniform  divisions  between  p = 0 and 
n = b/a.  The  mesh  point  positions  are  denoted  by  (C^^.Pj),  i = 1. 
j = l...m(i),  with  locating  the  mesh  point  at  the  discretized 

boundary  for  each  The  quantities  A?  and  Ap  denote  the  respective 

lengths  of  the  C and  p divisions,  with: 


• n; 


Ap  = (b/a)AC 


■ 

(4.1) 

. 


The  contact  pressure  distribution  is  modeled  by  assigning  a value  of  iji 
at  each  mesh  point; 

4>.  . , 1=1,... n ; j=l,...,m(i) 

» J 

with  (()  varying  linearly  in  C and  p between  adjacent  mesh  points. 

This  linear  variation  is  achieved  over  each  rectangular  element 
bounded  by  adjacent  horizontal  and  vertical  mesh  lines  by  dividing 
the  element  along  a diagonal  into  two  right  triangular  segments  over 
which  the  distribution  of  4>  is  planar.  This  is  shown  in  Figure  4-2 
for  a single  element,  and  Figure  4-3  shows  an  example  of  a complete 
pressure  distribution  comprised  of  such  elements.  Referring  to 
the  orientation  of  Figure  4-1,  the  dividing  diagonals  connect  the 
mesh  points  at  the  upper-left  and  lower-right  corners  of  each  element. 
Alternatively,  dividing  diagonals  connecting  mesh  points  at  the  lower- 
left  and  upper-right  corners  of  each  element  could  just  as  easily 
liave  been  used,  but  were  ruled  out  because  they  appeared  less  con- 
venient for  modeling  a contact  zone  boundary  of  the  general  shape 
shown  in  Figure  4-1. 
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Kig.  4-1  Discretization  of  Contact  /.one 


Fig.  4-2  Elt'mc'nt  i>f  I.inear  Fr 
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</>  (Cv^ 


V 


Fix-  0 i Kt ret  ized  Pressure  Distribution 
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II 
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The  discretized  form  of  Equation  3.11  may  now  be  written  for  the 
sum  of  the  normal  displacements  of  the  two  solids  at  each  contact 
zone  mesh  point  by  expressing  the  area  integral  over  the  contact 
zone  as  the  sum  of  the  area  integrals  over  the  component  elements 
of  linear  pressure  variation.  Figure  4-4  shows  the  (i,))  elemental 
rectangular  region  of  Integration.  The  linear  pressure  distributions 
over  the  right-triangular  subdivisions  labeled  "L"  and  "U"  are: 

. (C’.n’)  = AC  ,..-C+P(n'-n\)] 

i.j  t,j  1+1  j 


* * h,J+l  /« 


(4.2) 


where 


p H AC/All  = a/b 


(4.3) 


with  the  resultant  discretized  integral  equation  written  for  mesh 
points 


2 2 


n-1  m(l)-l  ^’j+1  C\+p(n'.^^-n') 

I ^ i f f AC,C 

i=l  j=l  n'.  ' 

J ^ I 


;’dn' 

(4.4) 


h j+l  ^ i+1 
+ 7 f 

n'.  C’i+P(n’j^l-n’) 


Q(^,,n,;f,’,n’)4’:'  , (f,’ .n’)d^'dn’ 

I J 1 • 1 


I«1 , . . . ,n  ; J»1 , . . . ,m(I) 


I 

I 


Fig.  4-4  F.Jemental  Region  of  Integration  for 
Discretized  Integral  Equation 
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Defining  the  quantities: 


■ L3 
.Y  J 


1 

AC 


j+1 

/ 


C'i+p(n'j^.i-n') 


^ j 


j+i-n  ; 

/ Q(Cj,nj;C',n')  jC[|  dc'dti’ 


^’l 


r u I 

a 


I U , 

Ly  j 


1 

AC 


j+1 


1+1 

/ 


I » 'J  > i * j 


Q(Cj,nj;C’,n’) 


1 

Cl 

Ln' 


(4.5) 


dC'dn' 


which  will  be  evaluated  by  numerical  Integration,  Equations  4.4  take 
the  form: 


r - 

r2  2 

- ""j  = 

n-l 

m(l)-l 

Z 

i;  {■)>. 

1=1 

j=i 

T • • - p(Yt  t • -“P'  -“t  t • C] 


u 


^l+l,j  f^I,J;l,j  ^'l“l,J;l,j  '^^'^'j+l“l,J;l,j  ^I,J;l,j^^ 

‘*’l,j+l  ^'l+l“l,J;l,j  ■ ®I,J;l,j  ^ 

'''  ‘^i  + l.i  + l j+l“l,J;l,j  " ^I,J;l,j^^' 


1=1 n ; J=l,...,m(I) 


I 
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By  adjusting  the  summation  indices  in  the  second,  third,  and  fourth 
lines  of  tin's  expression,  one  obtains: 


2 2 

r - 

1 

II 

n-l 

m(i)-l 

• ■ 

n m(  i 

-1) 

Z 

L 4 , 

,i=l 

+ 2 

Z 

i=l 

. I I , J; 1 , ] 

1 = 2 i 

=1 

n- 1 

m(i) 

I . -1 : 1 . J 

n m(i- 

1) 

r 

u 

i:  !>.  . 

j=3 

7.  1 

i=l 

i=2  j=2 

i . j I , i;  1 , ] 


i.j  I,.i;i,i 


(■4.6) 


where 


i = l , . . . ,n  ; ,1=1 , . . . ,m(I) 


,(1) 

,(2) 

,(3) 


(C',,1  + pn'.)  - PY^  . . 

i'hl  J I,,J,i,j  I,Jji,j 


^’i-l'^I.JU-l.j  ®I,,J;i-l,j  ^^'’’j+1  ^I,J;i-l,j  '^I.Jii-l,!^ 


'i',.,  51^  T.  ^ i 1 " T • • 1 " 0(n'i  ia 


(4.7; 


i-1  I,I;i,j-l 


i 1 1 pn'.)  . , . + . + 

I,J;i,j  i-1  j I,,l;i-l,j-l  I,.l;i-l,j-l 


PY 


U 

l,J;i-],,1-l 


It  is  assumed  that  the  boundary  contour  n = f(5)  is  monotonically 
decreasing  over  0 so  that  the  discretized  boundary  will  always 

satisfy: 

m(l)  = m(2)  = n 

m(l+l)  < m(i)  i=2,...,n-] 


Applying  this  condition  while  combining  the  four  double  summations 
in  Equation  4.6,  one  obtains: 
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2 2 

Cl  - enj 


n 

E 


m(i) 

E C 


i=l  j=l 


(4.9) 


I j J l,...,[n(I) 


with  the  influence  coefficients  C . given  by  Table  4-1  in  terms 

QS  ijJ » i>3 

of  the  quantities  C'  . defined  by  Equation  4.7,  The  subscripts 

J > i > J 

I,J;i,j  are  left  off  for  conciseness  in  Table  4-1. 


Equation  4.9  constitutes  N simultaneous  linear  equations  in  the  N 

discrete  pressures  (|)  . where 

i » J 


N = E m(i) 
i=l 


(4.10) 


For  given  values  of  the  elastic  properties  E^/E^,  E^/E^,  v^,  v^, 
and  v^,  geometric  ratios  a/h  and  a/b,  and  contact  zone  mesh  parameter 
n,  the  integrals  in  Equation  4.5  are  computed  and  tabulated  for  all 
mesh  coordinates: 


V I 1 » • • • * n y J— 1 y . ■ . y n 
i=l,...,n-l  ; j=l,...,n-l 


using  procedures  to  be  presented  in  Section  4.2.  Then,  from 
Equation  4.7  and  Table  4-1,  the  set  of  influence  coefficients  C,.  , . . 
may  be  obtained  for  contact  zones  of  any  discretized  boundary  shape 
satisfying  the  conditions  given  by  Equation  4.8. 


For  computational  purposes,  the  following  notation  is  adopted  for 
Equation  4.9: 


N 

- = E C ♦ 

<5  1 p;q  q 

p q=l  ^ 


with 


2 2 
= r - Cl  - 


p=l,...,N 


(4.11) 


(4.12) 
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TABLE  4-1 

INFLUENCE  COEFFICIENTS  FOR  EQUATION  4.9 


i 

j 

C 

Restrictions 

1 

1 

c^i) 

2, . , . ,n-l 

c(i)  ^ qO) 

n 

c(3> 

2, . . . ,n-l 

1 

c('>  + c(2) 

2, . . . 

c('>  + c^^)  ^ 

m(i+l)>2 

c^3)  + c^^) 

' 

m(i+l) 

c(2)  ^ (,(3)  ^ 

m(i+l) <m(i) 

m ( i+1 ) 

c(3)  -Hc('>' 

m(i+l)=m(i) 

1 

m(i+l)+l,.... 

C(2)  + 

in(i+l)  <m(  i)  -1 

in(i)-l 

m(i) 

c(^) 

m(i+l) <m(i) 

n 

1 

c(2) 

2, . . . ,m(n)-l 

c^2)  + 

m(n)>2 

m(n) 

c(^) 

(4.12) 
(Cont ' .d) 


C = C 

p;q 


1,1 


where  the  indices  p and  q replace,  respectively,  the  pairs  of 
indices  (I,J)  and  (i,j)  according  to: 


1=1; 

J=l, . . . ,n 

p = 

i-i 

[j  + E m(k) 

II 

M 

. ..,n  ; J=l, 

k=l 

1=1; 

a 

it 

q =< 

i-1 

i 

'j  + Z m(k) 

1=2, 

...,n  ; j=l. 

(4.13) 


k=l 


Since  any  discretized  boundary  allowed  by  Equation  4.8  contains 
the  vertices  (C  = 0,  q = b/a)  and  (C  = 1,  n = 0) , the  pressure 
is  implicitly  constrained  to  vanish  at  these  two  points.  The  criteria 
for  satisfactory  selection  of  the  boundary  shape  between  the  vertices 
is  that  the  discretized  pressure  distribution  obtained  as  a solution 
to  Equation  4.11  vanish  along  the  correct  boundary.  Clearly,  this 
boundary  selection  is  an  iterative  process. 

For  consistency  with  the  above-mentioned  constraint,  the  vertex 

pressures  (t  and  4 , ,,,  must  be  set  to  zero  in  each  of  Equations 

n n-m(n)+l 

4.11.  With  some  minor  rearrangement,  this  leaves  the  following  system 
of  N linear  equations  in  the  N-2  unknown  pressures  along  with  the 
unknown  parameters  T and  e: 

_ N _ _ 

6 = E p=l,...,N  (4.14) 

P p;q  q 

with 


p=l,...,N 


(4.15) 
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"r 

[ 

q = 

n 

q = 

N-m(n)+l 

1 

k 

all 

other  q 

r-‘ 

q = 

n 

c 

p;q 

(p=i , . , 

. .,N)  = 

U 

q = 

N-m(n)+l 

4 

p;q 

all 

other  q 

(4.15) 
(Cent  ' d . ) 


In  selecting  the  discretized  boundary  to  approximate  a continuous 
curve  as  shewn  in  Figure  4-1,  the  pressures  at  mesh  points  outside 
the  continuous  curve  are  set  to  zero  and  the  members  of  Equation  4.14 
prescribing  the  contact  surface  displacement  conditions  at  these 
points  are  effectively  discarded.  Letting  such  a point  be  denoted  by: 


(i  ,j  ) = (I  ,J  ) 
o o o o 


= 9, 


this  is  accomplished  by  setting: 

c . = r ’ “ 

^o’*^  l_0  all  other  q 


which  replaces  the  p^th  member  of  Equation  4.14  with  the  equation 

^ =0. 


The  resulting  system  of  Equation  4.14  is  solved  for  the  N i(>q's  using 
a standard  Gaussian  elimination  technique,  and  the  solution  is  in- 
spected for  approximate  satisfaction  of  the  condition  41  = 0 along 
the  selected  boundary  contour. 


The  d imens  ion  i e.ss  load  W,  defined  by  Equation  3.14,  is  obtained  from 
the  discretized  distribution  of  4>  by  forming  the  summation  over  the 
contact  zone  of  tiie  volumes  of  each  rectangular  pressure  element 
such  as  shown  in  Figure  4-2.  Letting  w^ 
an  element: 


'l.j 


denote  the  volume  of  such 


AO 


_ n-1 

W = A I Z 

i=l  j=l 


w 

i.  j 


where  it  is  clear  from  Figure  A-2  and  basic  concepts  of  solid 
geometry  that 


w.  . — T An  AC  (<j>.  .+2<t)  , .+2(|>.  1 . i) 

1.1  6 ' ^ '^i,j  ^i+l,j  ^i,j+l  ^i+l.j+l-^ 


(4 


Substituting  this  expression  for  w^  ^ into  Equation  A. 16,  one 
obtains  the  relation  which  may  be  used  directly  to  compute  W 
from  the  solution  to  Equation  A.IA: 


_ „ n-1  m(i-l)-l 

W = -IatiAC  z z 
1=1  j=l 


(((l.  . + 241  . . + 24  . ) 

1.1  1+1.1  1.1+1  1+1, J+1 


(4 


with  Equation  A, 13  used  to  convert  from  the  solution  notation  4 to 

_ q 

4.  ..  With  c,  r,  W,  and  the  4 distribution  available  from  this  dis- 
^ » J 

cretized  solution.  Equations  3.28  through  3.32  are  used  to  obtain 

the  quantities:  a /h,  a/a  , b/a  , d/d  , and  p/p°  . 

o o o o o 


A. 2 Evaluation  of  Influence  Coefficients  for  Discretized  Integral 
Equation 

A. 2.1  Formulation  in  Terms  of  Integral  Functions 

The  integrals  appearing  in  Equation  A. 5 of  the  function. 


Q(Cj.nj;C',n')  = K(Cj-C' .hj-n’)  + K(Cj+C' .Oj-n') 
+ K(Cj+C'.nj-n')  + K(Cj+C'.nj+ri') 


defined  by  Equation  3.12,  must  be  evaluated  numerically  in  order  to 
determine  the  influence  coefficients  C from  Equation  A. 7 and 

Table  A-1.  Equation  A. 5 can  be  broken  down  into  the  integral  functions. 


1-' 

n 

C _ ' 

(C.n)  = / 

/ K(u,v) 

dudv 


(4 


.16) 


.17) 


.2.1) 


I 


A1 


[ H*1 

1 

n 

1 

> 

a 

\ 

I* 

* (C,n,+p)  = / 

1 

S' 

c 

< 

' dudv 

(4.2.1) 

/ 

i 

i 

i] 

1 

(Cont ' d . ) 

which  are  naturally  amenable  to  systematic  numerical  computation  and 
tabulation.  To  express  Equation  4.5  in  terms  of  these  integrals, 
first  substitute  the  above  expression  for  Q into  the  first  three  of 
Equation  4.5  and  change  the  variables  of  integration  in  each  term 
such  that  they  are  the  same  as  the  arguments  of  K to  obtain: 


L-l 


n,-n 


t.Y 


J I.J:i,j 


AC 


n .-n  . 
•J  J 


[C,-C'  .+P (n.-n -pv 


J ' j+1  i 'J  ' j+1' 


f 

V^’i 


-1' 

K(u,v)  • dudv 


v"'j+i 


AC 


n.-h  . 
J J 


n.+  n . , , 

J J+1 


* V”'- 

J J 


^ ‘5  n/r,'. 
J J 


v^'i 


' ] 

K(u,v)  • 


;-n  \ 
'j> 


[Cj-C'..-I 


f K(u,v)->u-c 


V^’i 


(4.2.2) 

dudv 


+p(n,+n' . , , ) ]-pv 

J J+i  i i J j+1 

/ K( 


u,v) • ju 


-'l, 

iy-Pj, 


dudv 


It  should  be  clear  that  this  expression  can  be  applied  to  the  latter 
three  of  Equation  4.5  by  changing  the  sign  of  each  integral  and 
replacing  C'^  with  the  lower  limit  of  the  u integration  for 

each  term.  Now  separate  Equation  4.2.2  into  integrals  over  rectangular 
regions  denoted  as: 


fAl 


[ 


V"'j+i  ^r^'i  _ 


J 


K(u,v)  <Ct-u) 


V^'j+1 


n -n ' 

J j 


- r‘ 

K (u,v)  tC^-u 


dudv 


(4.2.3) 


dudv 
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v^’i 

- / 

f 

K (u.v)- 

o 

V^'j+1 

v^'i 

- / 

/ 

K (u.v) 

n ,+n ' . 

J 3 

0 

integrals  over 

trapezoidal  regions 

denoted 

fA*! 

h T“h  ' -.Li 

J j+1 

[Ci-C’i+P(nj-n'.^^ 

- 

/ 

f 

I,J;l,j 

v^'j 

o 

-1  1 


dudv 


(4.2.3) 
(Cont 'd . ) 


K(u,v)  • dudv 


V'^'j+l  [Ci+C'.-P(nj-n’.^^)]+pv 


+ / 


n.-P  • 
J 3 


f 

O 


K(u, 


V^’j+i  tv^’rp(v^’j+i)J+pv 

+ f , f 

1 .n' 

j""  J 


K(u, 


(4.2.4) 


V^’j+1 


n ,+ri ' . 
j 3 


[Cj+C ' ^+P  (Pj+n  V^^)  ]-pv 

/ K (u,v)  • 

o 


to  obtain  the  following  convenient  representation  of  all  six  of 
Equation  4.5: 


' r L-i 

\ ® 

A* 

A* 

) 1 

1 

1 

<6 

1 1 ' 

° AC  ' ' 

11 

+ TT—  • 

AC 

1 

1 ^ 

fc 

G* 

- 

(4.2.5) 


U) 


1 

AC 


1 

AC 


i.J;i+i.j 


A* 

B*) 

G* 


Equation  4.2.3  for  A,  B,  and  G can  be  written  in  terms  of  the  integral 
functions  H,  I,  and  J to  yield: 


' n J 


+ - -f  CjH-l 


- U"-'  v’Vi'  -"k"-']  v^'j* 


(4.2.6) 


l-CjH  (Cj-C'i,n/n’.^,)  + l-CjH  (Cj-C',.n^+n’.) 


njH-j 


i-?iH  (Cj-bC'i.n/n'.^l)  + .C,+C',.n,+n’.) 


_J-njH 


. J-njH  > 


Similarly  A*,  B*,  and  G*  may  be  expressed  in  terms  of  H*,  I*,  and  J*: 


i,d; i,  j 


j*-n  jH* 


(4.2.7) 


H*  ' H*  1 
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with 


.(1) 


.(2) 


= - C.  + P(nj-n'.^^) 


C-  = C,  + C-,  - P(nj-n’ 


= Cj  - - P(nj+n'.^^) 


£<“>  + »(V"'j+l> 


(4.2.8) 


Equation  4.2.5  together  with  Equations  4.2.6  and  4.2.7  give  the 
desired  expression  of  the  quantities  a^,  3^,  a^,  3^,  and 

in  terms  of  H,  I,  J,  H*,  I*,  and  J*. 


Recalling  from  Equation  3.7  that  K(C,n)  = aK(aC.an)  with  a = a/h 
and  K(u,v)  the  kernel  function  given  by  Equation  3.5,  it  is  seen  that 
H,  I,  J,  H*,  I*,  and  J*  are  obtainable  from  computations  of  the 
basic  integral  functions: 


ri 

t 

1 

(s,t)  = f 

f K(u,v) 

h‘ 

1 

Lv 

(4.2.9) 


r H*' 

1 ) 

\ 1 

t s+pv 

I*' 

(s,t,+p)  = / / K(u,v) 

1 4 

o o 

^ V ] 

dudv 


with 


H(£;,ri)  = (1/a)  H (a^.an) 

I(Cin)  = (1/a^)  I (aC,an) 

J(Cin)  = (1/a^)  J (aC.an) 
H*(f,,ri,+p)  = (1/a)  H*  (ar,,ar),+p) 
I*(C,n,+p)  = (1/a^)  I*  (af,,(tn,+p) 
“ (1/a^)  J*  (a.f;,(«n,+p) 


(4.2.10) 


I 


From  Equation  3.5,  it  is  clear  that: 


K(u,v)  = K(-u,v)  = K(u,-v)  = K(-u,-v) 
and  as  a consequence  of  this  symmetry: 

H(s,t)  = Sgn(st)  H (|sl , |t|) 

I(s, t)  = Sgn(t)  I ( I s! , I t 1 ) 

J(s,t)  = Sgn(s)  J (Isj.lt!)  (4.2.11) 

H*(s,t,+f;)  = Sgn(st)  H*  ([s|,!t|,  +p  Sgn(st)] 

I*(s,t,+p)  = Sgn(t)  I*  [is|,|t|,  +p  Sgn(st)] 

.I*(s,t,+p)  = Sgn(s)  J*  [|sl,|t|,  +p  Sgn(st)l 

where 

i-l  X < 0 

(4.2.12) 

+1  X ^ 0 

so  that  the  basic  integral  functions  defined  by  Equation  4.2.9  need 
only  be  computed  for  positive  values  of  s and  t.  By  combining 
Equations  4.2.10  and  4.2.11,  one  obtains  the  relations: 

H(C,n)  = (l/ot)  Sgn  (Cn)  H (a I f,  I ,a  I n I ) 

i(?,n)  = (l/ri^)  Sgn  (n)  J (fi|  cl  ,^i|  nl ) 

J(c,n)  = (1/a^)  Sgn  (c)  J (a|c|  ,ni|n!)  (4.2.13) 

H*(C.n,+p)  = (l/ot)  Sgn  (Cn)  H*  [m  | C | ,<i  1 n | ,+  pSgn  (Cn)l 

I*(C,n,+p)  = (l/oi^)  Sgn  (n)  I*  fa|  C I ,a)  n|  ,+p  Sgn  (Cn)] 

; J*(C,n,+p)  = (1/ot^)  Sgn  (C)  J*  [a  | C 1 .'’t  | n i .+pSgn  (Cn)  1 

which  will  be  used  to  obtain  the  values  of  H,  I,  J,  H*,  I*,  and  J* 
required  in  Equations  4.2.6  and  4.2.7  from  computed  values  of 
H,  I,  J,  H*,  I*,  and  J*. 

' Equations  4.2.6,  4,2.7,  and  4.2.8,  together  with  Equation  4.2.13,  must 

^ now  be  closely  examined  to  determine  all  arguments  s,t  for  which  H,  I, 

J,  H*,  1*,  and  J*  must  be  evaluated  in  order  to  provide  the  values  of 
H,  I,  J,  H*,  I*,  and  J*  at  all  arguments  called  for  in  Equations 

I 
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4.2.6  and  4.2.7.  First,  note  that  the  quantities  (A.B.G) 
given  by  Equation  4.2.6  and  (A*,B*,G*) 


I. J;  i,  i 
given  by  Equation  4.2.7 


must  be  obtained  for  all  the  Integer  subscripts  1=1,..., n;  ./=],... n; 

i=l,...,n-l;  i=l,...,n-l,  in  order  to  obtain  all  of  the  Influence 

coefficients  C , from  Equation  4.2.5,  Equation  4.7,  and  Table  4-1 

^ ^ ^ * J 

Thus,  the  mesh  coordinates  6', ,n'.,  will  take  on  their  full 

I J 1 j 

range  of  values: 


(I-l)AC 

1=1,.. 

. .n 

(J-l)An 

J=l,.. 

. ,n 

(i-l)AC 

1=1,.. 

(j-l)An 

j=l,.. 

. ,n 

in  forming  the  arguments  for  which  H,  I,  J,  H*,  I*,  and  J*  will  be 
needed  in  Equations  4.2.6  and  4.2.7.  Clearly,  for  any  values  of 
[H,  I,  J)  (C/i)  and  fH*,  I*,  J*]  (C,n,  +p)  required  in  Equations  4.2.6 
and  4.2.7,  the  argument  C will  be  a positive  or  negative  Integer 
multiple  of  AC  and  the  argument  q will  be  a positive  or  negative 
integer  multiple  of  An.  For  the  corresponding  values  of  [H,  I,  J] 
(s,t)  and  [H*,  I*,  J*]  (s,t,  +p)  called  for  in  Equation  4.2.13,  the 
arguments  s and  t will  be  positive  Integer  multiples  of  aAC  and 
aAn,  respectively.  Computations  of  [H,  I,  J]  (s,t)  and  [H*,  I*,  J*] 
(s,t,  +p)  are  thus  required  over  particular  ranges  of  s at  intervals 
of  aAC  and  t at  intervals  of  aAn.  By  noting  the  extreme  values  of 
the  arguments  of  H,  I,  J,  H*,  I*,  and  J*  in  Equations  4.2.6  and  4.2.7 
over  the  full  range  of  values  of  the  mesh  coordinates,  and  using 
Equation  4.2.13  to  determine  the  corresponding  extreme  values  of  the 
arguments  of  H,  I,  J,  H*,  I*,  and  J*,  one  can  establish  the  following 
ranges  on  s and  t for  computation  at  respective  intervals  of  aAC 
and  aAn: 


[H,I,J1  (s,t) 

(.s.t,+p) 


r 0 ^ s ^ 2(n-l)aAC 
/.  0 < t ^ 2(n-l)aAn 

0 ■ s ^ 2(n-l)aAC-Pt 

.0  ■ t £ (n-1 )aAn 


(4.2.14) 
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(s.t.-p) 


[0  ^ s ^ 4(n-l  )fjAC 
0 ^ t < (n-l)aAr) 


-pt  - (n-l)aA/^  i “ H 4(n-l)aAA, 
■ nriAn  < t < 2(n-])''tAri 


(4.2.14) 
(Cont ' cl . ) 


Ely  utilizing  the  computational  methods  derived  in  Section  4.2.2, 

the  desired  arrays  of  values  of  (s,t)  and  [11*, I*, J*]  (s,t,+p) 

are  obtained  over  these  ranges.  These  numerical  arrays,  together 

with  Equations  4.2.13,  4.2.6,  4.2.7,  4.2.5,  and  4.7,  and  Table  4-1, 

provide  for  the  complete  determination  of  the  array  of  Influence 

coefficients  C . in  the  discretized  Integral  equation  given 

1 » 1 1 i > 1 

by  Equation  4.9. 


4.2.2  Evaluation 


Func  t ions 


Substituting  from  Equation  3.8  for  the  approximate  form  of  K(u,v) 
in  the  Integral  functions  defined  by  Equation  4.2.9  and  separately 
defining  the  Integrals  of  the  Hankel  inversion  integral  and  inverse 
square  root  parts  of  K(u,v),  one  obtains: 


I \ (s.t)  yj  I 


I * ''  ( s , t , +p ) 


(s,t)  + (1+6) 


(s,t,+p)  + (1+6) 


(s,  t) 


(s,t,+p) 


(4.2.15) 


1 ( (s,t) 


(.s,  t ,+c> ) 


J'  [g2(“)'g^(M)-l]  / -^0^“  ) 1 u ^ dudvdu) 

o o o \ 

v 


(4.2.16) 


lgT(  • )-g,  ('■^)-l  1 / J (cii'^  +vS  1 u ( dudvdw 


I 


48 


and , 


t 

® 2 

lo 

, 

(s,t)  = 

f 

f (u  +v 

J 1 

o 

o 

. 0, 

\ 

Ch*] 

t s+pv 

I*' 

• (s,t,+p) 

= 

/ / 

/ J* 
o- 

0 o 

2, -1/2 


dudv 


, 2.  2, -1/2 

(u  +v  ) 


’1 


dudv 


(4.2.17) 


where  in  Equation  4.2.16,  the  natural  order  of  integration  has  been 
switched  by  moving  the  w integration  to  the  outside. 

The  u-v  double  integration  in  Equation  4.2.17  is  performed  analytically 
by  following  the  procedure  described  in  Appendix  I.  Defining  the 
quantities: 


a = 


-/l  + 


R- 


/2  . 2 

/s  + t 

s + pt 

/+2~~~2 
✓ s—  + t 


the  resulting  closed-form  expressions  in  s and  t take  the  form: 


^ , R"  + t . ^ , R°  + s 

H^(s,t)  = slog  ^ +t  log  


1 9 R®  + ^ 1 

I (s,t)  = -J  s log  + tCR'-t) 

J^<s,t)  ■=  ^ s(R“-s)  + ^ log  ^ ^ 

. ? log 

O - O T 

OS-  + pR" 

. _ / +.  . R""  + I s — I 

+ Sgn(s-)t  log  


(4.2.18) 


(4.2.19) 


I 

I 
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I*  (s,t,+p) 


J*  (s,t,+p) 


i (t) ' F 


+p(R-  - s) 


I t(R^  - t) 


i (5) 


+ i log 

a ^ + - + 

as-  + R- 


+ P log 

s a 

as 


(ot  + R-)  ' 

i » Rt  J 


] 


(4.2.19) 
(Cont ' d . ) 


^+.1  2 ^ R-  + 

+ Sgn(s-)  y t log  — 


The  u-v-w  triple  integration  must  be  carried  out  numerically  in 
evaluating  the  functions  defined  by  Equation  4.2.16.  It  is  seen  that 
each  s,t  argument  pair  within  the  ranges  given  by  Equation  4.2.14 
defines  a region  of  u-v  integration  of  one  of  six  types  shown  in 
Figure  4-5.  Decomposing  these  regions  in  the  manner  shown,  a 
numerical  scheme  is  derived  for  evaluating  the  functions  at  Intervals 
of  aAC  in  s and  aAn  in  t.  For  Region  I.  this  decomposition  allows 
[Hj^,  I^,  (s,t)  to  be  expressed  in  the  form: 


i:;! 

Ki 

‘ ^1/  (s,t-aAn) 


[g2(ui)-g^(“;-il 


t s 

/ / 

t-aAn  s-aAC 


J (w, 
o 


+v' 


L-\ 


dudvdu) 


(4.2.20) 


I,)  (s-aAC,t)  - 


I 

J 


(s-aAC , t-aAn) 


Considering  Regions  II  through  V,  [H*.  I*,  J*]  (s,t,-p)  can  be 
expressed  in  a similar  manner  for  the  various  ranges  of  s: 


s = pt 


fH*! 


C'tj 


y (pt,t,-p)  = 


(4.2.21) 
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»1 

+ }1* 


(pt,  t-aAn,-p) 
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(Cont'd.) 
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Finally,  for  Region  VI: 


/_J* 


(s,t,+p) 


o t s+pv 

/ [g  f ^ 

o 2 ^ t-aAn  s-aAC+pv 


m* 


(s, t-aAn, +p)  + 


H*' 


L 
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( s-aAC , t-aAn ,+p) 


In  addition,  note  the  following  from  inspection  of  Equation  4.2.16: 


fH*')  ' 

ol 

\ 1 

\ 1 

(s,0)  = I^ 

, (0,t)  = jl*i  (s,0,+p)  = 

r 

/ -^1 

, /j*( 

0 

t 1 . 

( 1 - 

L Ij 

(4.2.26) 


T-H*  1 

j I*  > (0,t,-p) 

L-^\\ 


(4.2.27) 


Equations  4.2.20  through  4.2.27  are  applied  directly  in  evaluating  the 
functions  at  successive  increments  of  aA^  in  s and  aAn  in  t between 
the  limiting  values  given  in  Equation  4.2.14.  Each  function  evaluation 
entails  a numerical  quadrature  over  the  shaded  segment  of  one  of  the 
regions  shown  in  Figure  4-5,  with  quadrature  over  the  remainder  of  the 
region  having  occurred  for  function  evaluations  at  preceding  values 
of  s and  t. 


Since  the  Integrands  requiring  numerical  quadrature  are  smoothly 
varying  continuous  functions  of  lO,  u,  and  v,  standard  single-variable 
Gaussian  quadratures  are  used  to  integrate  in  each  variable  in  the 
order  shown.  To  avoid  any  numerical  difficulties  arising  from  the 
oscillation  of  J , the  u)  region  (0  to  u'  ) is  divided  into  uniform 


25) 


intervals,  roughly  equal  in  length  to  a period  of  J (Rn)  where  R 

f?.  T ° 

is  an  estimated  mean  value  of  Wu  +v  over  the  u-v  region.  A 
ten-point  quadrature  in  oj  is  then  used  to  integrate  over  each  interval  . 

At  each  quadrature  abscissa  in  u,  [ g^ (to)-g^ (w) ] is  obtained  from  Equations 
2.14  through  2.16,  and  the  u-v  integration  is  performed  with  a quadra- 
ture in  u providing  the  value  of  the  integrand  at  each  abscissa  of 
a quadrature  in  v.  Polynomial  approximations  given  in  Reference  [ 9 ] 
are  used  to  compute  to  eight-place  accuracy  in  the  u quadrature. 

It  is  seen  from  Equations  4.2.20  through  4.2.25  that  the  u-v  limits 
of  integration  span  regions  of  at  most  aA^  and  aAn,  respectively. 

For  all  values  of  a,  AC,  and  An  considered  in  this  study,  and,  with 
03^  = 10,  the  argument  of  always  varied  less  than  a period  between 
the  u-v  limits.  Consequently,  convergence  to  seven  places  was  attain- 
able by  dividing  the  w region  as  described  above  and  using  5-point 
quadratures  over  the  full  u and  v regions.  Treating  the  oscillation 
of  as  approximately  sinusoidal,  the  number  of  ou  divisions,  n^ , was 
determined  from  the  formula: 

w R 
o 

''d  " 2 TT 

with  the  constraint,  n^  ^ 6.  Tables  of  abscissas  and  weight 
factors  used  for  the  5-  and  10-point  quadratures  are  given  in 
Reference  I9]. 

It  should  be  noted  that  for  fixed  AC  and  An,  with  AC  = l/(n-l)  and 
An  = AC/p  from  Equation  4.1,  the  magnitude  of  each  value  of  s and  t 
for  which  the  integral  functions  are  to  be  evaluated  is  proportional 
to  the  ratio  a = a/h.  Thus,  each  R and  consequently  each  n^  is 
also  proportional  to  a;  and,  since  the  number  of  abscissa  points  in 
, w between  the  limits  0 and  is  10  n^,  the  computing  time  required 

to  evaluate  each  triple  integral  is  similarly  proportional  to  a. 

T This  is  an  important  consideration  in  planning  a parametric  study  of 

' tlie  effect  of  a on  the  behavior  of  the  contact  zone. 


I 
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Once  all  the  values  of  H^,  I^,  J^,  H*,  I*,  and  J*  have  been  computed 

and  tabulated  using  the  above  methods,  they  are  combined  according 

to  Equation  4.2.15  with  the  corresponding  values  of  H , I , J , H*, 

o o o o 

I*,  and  J*  computed  directly  from  Equations  4.2.19  and  4.2.18,  to 
give  the  desired  array  of  values  of  the  integral  functions  H,  I,  J, 
H*,  I*,  and  J*  needed  to  form  the  quantities  considered  in  Section 
4.2.1  above . 


At  this  point,  it  is  worth  noting  that  the  elastic  properties  of  the 
homogeneous  solid,  E^/E^  and  v^,  enter  the  numerically  formulated 
contact  problem  only  in  Equation  4.2.15  through  the  parameter  B de- 
fined by  Equation  3.6.  The  elastic  properties  of  the  layered  solid, 
^1^^2’  '^l’  ^2'  the  functions  H^,  I^,  , H*,  I*,  and 

J*  through  g2(<^)-g^(w)  given  by  Equations  2.14  through  2.16. 


4 . 3 Computer  Programs 


Two  FORTRAN  programs,  HIJPRG  and  PS0LV,  implement  the  computational 
scheme  derived  in  Sections  4.1  and  4.2.  Both  programs  are  written 
for  a CDC  6600  computing  system  with  the  NOS  Operating  System.  Input 
and  output  descriptions  and  listings  are  given  in  Appendices  III 
through  VI. 


In  HIJPRG,  the  arrays  of  integral  functions  , 1^^,  , H*,  I*,  and  J*, 

computed  by  numerical  quadrature  are  obtained  for  prescribed  values  of 
E^/E2,  Vj^,  V2,  a,  p,  and  n as  described  in  Section  4.2.2.  In  PS0LV, 
the  Influence  coefficients  ^ ^ . are  generated  from  these  arrays 
through  Equations  4.2.19,  4.2.15,  4.2.13,  4.2.6,  4.2.7,  4.2.5,  and  4.7 
and  Table  4-1.  With  Equation  4.9  expressed  in  the  notation  of 
Equation  4.14,  a standard  matrix  inversion  subroutine  yields  the 


dimensionless  contact  pressure  distribution  (fi  .,  approach  P,  curva- 

^ * J 

ture  ratio  e,  and,  through  Equation  4.17,  load  W.  Finally,  the 

quantities  a /h,  a/a  , b/a  , d/d  , and  p°/p°  are  obtained  from  Equations 
0 0 00  0 

3.28  through  3.32.  The  analysis  contained  in  PS0LV  pertains  to  pre- 


scribed values  of  E^/E^,  v^,  boundary  shape  m(l)  1=1,..., n,  and  boundary 
points  (i,j)  of  zero  pressure. 
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PART  5 


NUMERICAL  RESULTS 


5.1  Test  for  Numerical  Accuracy 

With  the  integral  functions  H*,  I*,  and  J*  computed  to 

seven-decimal-place  accuracy,  the  numerical  formulation  is  grounded 
on  the  essentially  exact  computation  of  surface  d isplacements  due  to  the 
discretized  pressure  distribution.  The  error  in  the  solution  is  due 
to  the  truncation  error  in  approximating  the  true  continuous  pressure 
distribution,  an  error  that  in  principle  can  be  made  as  small  as  desired 
by  using  a fine  enough  mesh  to  construct  the  approximate  distribution. 

To  keep  computer  time  and  storage  at  moderate  levels  in  the  present 
study,  the  mesh  was  limited  to  ten  divisions  along  the  major  and  minor 
halfwidths  (n=ll).  The  truncation  error  associated  with  this  mesh  will 
now  be  examined  by  comparing  numerical  predictions  for  the  limiting 
cases  of  elliptical  Hertz  contact  of  homogeneous  solids  and  axisymmetric 
contact  of  layered  solids  with  available  exact  solutions.  These  compari- 
sons serve  as  an  overall  check  on  the  validity  of  the  numerical  formula- 
tion. 

As  discussed  in  Section  3.2,  the  formulation  for  layered  contact 
reduces  to  the  Hertz  problem  for  homogeneous  solids  when  one  prescribes 
E^/E2  = 1 and  = v^.  This  reduction  occurs  in  the  numerical  formula- 
tion through  Equations  4.2.15  and  4.2.16,  where  the  integral  functions 
Hj^ , I^,  Jj^,  H*,  1*,  and  J*  vanish  identically  for  all  arguments.  With 
n=ll,  the  elliptical  boundary  for  any  ellipticity  ratio  p is  prescribed 
by: 


i 

1 2 

3 4 5 6 

7 

8 9 10 

11 

m(i) 

11  11 

11  11  11  10 

10 

9 8 7 

5 

Tliis  is  the  boundary 

exhibited  by  the 

example 

contact  zone 

mesh  shown 

in  I’iguri 

' 4-1.  TIk' 

pressure  is  set  to 

zero  at 

the  following  boundary 

mesli  points 

wliieh  fall  outside 

the  true  elliptical 

contour : 

= (2,11) 

(3,11)  (4,11) 

(5,11) 

(6,10) 

(7,10) 

(7,9)  (8,9) 

(9,8) 

(9,7) 

(10,7) 

(10,6)  (11,5) 

(11,4) 

(11,3) 

(11,2) 
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With  E^/E^  = 0 completing  the  definition  of  this  approximate  licrtz 
problem,  numerical  solutions  were  obtained  at  ellipticity  ratios  of 
1,  2,  5,  and  10.  Table  5-1  shows  the  resulting  values  of  c,  F,  ((>°, 

W,  s/a^,  p°/p°,  and  along  with  the  exact  values  given  by  Equations 

3.18  through  3.23  and  3.33.  Note  the  improved  agreement  in  a/a  and 

d/d  relative  to  the  corresponding  quantities  W and  F.  This  is  ex- 

° — 1/3 

plained  in  part  by  the  proportionality  of  a/a^  to  W tiirough 

Equation  3.30. 


A second  comparison  is  provided  by  an  essentially  exact  computation  of 

surface  displacements  for  a layered  solid  due  to  an  axisymmetric 

ellipsoidal  pressure  distribution.  In  terms  of  the  dimensionless 
% 2 

variables  u = u R /a  for  surface  displacement  and  X = r/a  = 
o o X 


for  radius.  Equation  II. 1 of  Appendix  II  takes  the  form: 


u (X)  = — z — r — / [g.,(w)  - g,  (w)  - 1]J  (aXw)  (sinaw  - aucosaw) ^ (5. 

o ALtt  24  o 

o (aw) 

+ 2 - X^J  , X < 1 

for  the  surface  displacement  due  to  pressure  distribution: 

X £ 1 
X > 1 


For  a given  value  of  X,  the  Integral  in  Equation  5.1  is  evaluated  to 
7-place  accuracy  by  dividing  the  truncated  region  of  integration  (0  to  10) 
into  a suitable  number  of  intervals  as  discussed  in  Section  4.2.2  and 
using  a 10-point  Gaussian  quadrature  to  integrate  numerically  over 
each  Interval.  This  exact  computation  of  a^(X)  is  compared  with  the 
result  obtained  by  computing  the  right-hand  side  of  Equation  4.9  for: 


+ 1 1 


with  influence  coefficients  obtained  for  p = 1,  8=0,  and  the  m(i)'s 
given  above  for  a general  elliptical  boundary.  Computations  were  per- 
formed for  the  case  of  an  incompressible  layer  (v^^  = 0.5)  on  a rigid 
substrate  (E^^/E^  = 0)  with  a = L,  2,  and  4.  (fi®  was  assigned  a value 
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TABLE  5-1 

COMPARISON  BETWEEN  NUMERICAL  AND  EXACT 
SOLUTIONS  FOR  HERTZ  CONTACT 


p 

1 

2 

5 

10 

e 

1.000 

2.855 

11.95 

37.05 

Computed 

1.000 

2.843 

11.83 

36.54 

1 Exact 

0.0 

+0.4 

+1.0 

+1.4 

% Error 

r 

1.96 

1.68 

1.46 

1.36 

Computed 

i 

2.00 

1.71 

1.47 

1.37 

Exact 

-2.0 

-1.8 

-0.7 

-0.7 

% Error 

0.402 

0.502 

0.777 

1.180 

Computed 

0.405 

0.505 

0.777 

1.176 

Exact 

-0.7 

-0.6 

0.0 

+0.3 

% Error 

w i 

0.825 

0.515 

0.320 

0.243 

Computed 

i 

0.849 

0.529 

0.326 

0.246 

Exact 

1 

1 

-2.8 

-2.6 

-1.8 

-1.2 

% Error 

a/a 

o 

1.010 

1.181 

1.385 

1.517 

Computed 

1.000 

1.171 

1.376 

1.511 

Exact 

+1.0 

+0.9 

+0.7 

+0.4 

% Error 

pvp; 

1.002 

1.464 

2.656 

4.42 

Computed 

1.000 

1.459 

2.640 

4.38 

Exact 

+0.2 

+0.3 

+0.6 

+0.9 

% Error 

d/d 

0 

1.000 

1.173 

1.396 

1.560 

Computed 

1.000 

1.173 

1.395 

1.557 

Exact 

0.0 

0.0 

+0.1 

+0.2 

% Error 
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of  4/tt  to  normalize  the  coefficient  of  Equation  5.1.  Table  5-2  shows 
the  resulting  exact  and  numerically  computed  values  of  for  \ 

between  0 and  1 in  increments  of  0.2.  As  expected  from  the  nature  of 
the  discretized  pressure  distribution,  agreement  is  good  near  the 
center  (A  = 0)  and  relatively  poor  at  the,  edge  (A  = 1)  where  the 
gradients  of  the  actual  pressure  distribution  are  high  and  not 
adequately  represented  by  the  discretized  model.  The  difficulty  at  the 
edge  is  obvious  from  an  inspection  of  Figure  4-3  which  illustrates  the 
discretized  model  of  an  ellipsoidal  pressure  distribution  for  n=]l. 

The  final  comparison  is  provided  by  the  solutions  reported  by  Chen 

and  Engel  [4]  for  axisymmetric  layered  contact.  As  discussed  in 

Part  1,  these  are  approximate  solutions  to  a one-dimensional  integral 

equation.  Nevertheless,  the  accuracy  level  verified  by  the  authors  is 

high  enough  that  the  solutions  can  be  assumed  exact  for  purposes  of 

this  comparison.  The  reported  solutions  consist  of  tabulated  values  of 

a dimensionless  load  p*  = 3ttW/8  and  approach  6*  = 1/2  over  ranges  of 

a and  with  ~ ® '^l  ” ^^3  ” Values  of  W and  F 

computed  by  the  present  method  for  a = 1,  5,  = 0.1,  10  are  shown 

in  Table  5-3  along  with  the  corresponding  reported  values.  Also  shown 

in  this  table  are  the  values  of  a/a  and  d/d  obtained  from  W and  F 

o o 

through  Equations  3.30  and  3.28.  Relative  error  magnitudes  are  notably 

similar  to  those  shown  in  Table  5-1  for  contact  of  homogeneous  solids.  ; 

Again  note  the  small  errors  in  a/a  and  d/d  relative  to  those  in  W 

o o 

and  F . ; 

I 


The  comparisons  with  the  Hertz  solution  and  the  axisymmetric  analysis 
of  Chen  and  Engel  show  an  error  of  roughly  1 percent  in  the  solution 
quantities  of  primary  physical  Importance:  curvature  ratio  c,  major 
halfwidth  approach  d/d^,  and  central  pressure  p°/p°.  The  edge  dis- 

crepancies in  the  comparison  between  exact  and  computed  displacement 
profiles  are  representative  of  local  errors  that  can  be  expected  in 
computed  pressures  near  the  edge  of  tlie  contact  zone.  Tlius,  the 
present  formulation  with  n=ll,  while  providing  a fairly  accurate  computa- 
tion of  e,  a/a^,  d/d^,  and  p'^/p°,  does  not  given  an  accurate  prediction 
of  either  the  pressure  profile  shape  near  the  boundary  or  the  shape  of 
the  boundary  contour.  A finer  mesh  would  be  required  to  accurately 
determine  these  shapes. 
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TABLF.  5-2 

COMPARISON  .f5ETWEEN_NUMERI_CAL  .ANAA^’^ 
spji^JJA^A_APA.AiA/ACE_  piPA'APAMEM'rp  due  to 

AXISYMMETRTC  ELLIPSOIOAI.  PRESSURE  DISTRIBUTION 


1 

* ct 

mim 

1 

2 

4 

c 

II 

r< 

c 

<r  - 

0.755 

0.264 

0.0579 

Computed 

2.00 

0.755 

0.263 

0.0574 

Exact  i 

1 

i 

0.0 

+0.4 

+0.9 

X Error 

‘ >=0.2 

0.732 

0.261 

0.0583 

Computed 

1.96 

0.732 

0.260 

0.0578 

Exact 

1 ” 

0.0 

+0.4 

+0.9 

% Error 

1 A=0.4 

1 

1 - 

0.659 

0.247 

0.0593 

Computed 

1 

1 1.84 

0.660 

0.246 

0.0587 

Exact 

- 

-0.2 

0.4 

+1 . 0 

% Error  | 

>=0.b 

1 

- 

0.0591 

Computed  j 

1 

1.64 

0.0583 

Exact  1 

- 

+ 1.4 

% Error  i 

>=0.8 

■1 

0.342 

0.130 

0.0453 

1 

Computed  1 

0.348 

0.133 

0.0456 

Exact 

\ 

■■ 

-1.7 

-2.3 

-0.7 

% Error 

> = 1.0 

1 

0.077 

-0.030 

■ 

-0.0356 

Computed  | 

1.00 

0.090 

-0.020 

-0.0286 

Exact 

- 

-14 

+50 

+24.5 

% Error  1 

Elastic  Properties:  = 0 

Vj^  = 0.5 


Pressure  Distribution:  4i 


= [(4/71^) 


\ < 1 
> > 1 


5.2  Numerical  Solutions 


Numerical  solutions  over  ranges  of  a and  p have  been  obtained  for  three 
sets  of  elastic  property  ratios  in  order  to  determine  the  influence  of 
layer  thickness  and  curvature  ratio  on  the  major  and  minor  halfwidths, 
approach,  and  central  pressure  for  a given  load.  The  elastic  property 
ratios  pertain  to  (i)  contact  between  a rigid  solid  and  an  incompressible 
layer  backed  by  a rigid  substrate,  (ii)  contact  between  a rigid  solid 
and  a high  modulus  layer  backed  by  a low  modulus  substrate,  and  (iii)  con- 
tact between  a low  modulus  solid  and  a high  modulus  layer  backed  by  a 
low  modulus  substrate.  The  following  values  represent  these  configurations 


Ei/E^ 

E1/E3 

"1 

"2 

"3 

(1) 

0 

0 

0.5 

- 

- 

(ii) 

10 

0 

0.25 

0.25 

- 

(iii) 

10 

10 

0.25 

0.25 

0.25 

Computations  were  performed  with  values  of  1,  2,  3,  4,  and  5 assigned 
to  a and  with  values  assigned  to  (i  in  increments  of  0.5  from  1.0  to 
the  first  value  resulting  in  a computed  curvature  ratio  e of  10  or  above. 
In  all  cases,  a discretized  elliptical  boundary,  specified  according 
to  the  values  of  m(l)  and  points  of  zero  pressure  given  in  Section  5-1 
for  n=ll,  resulted  in  a qualitatively  well  behaved  pressure  distribution 
in  the  sense  that  pressures  fell  to  zero  at  the  boundary  without  any 
unexpected  Irregularities.  In  other  words,  any  departure  of  the  boundary 
from  the  ellipse  of  the  Hertz  solution  could  not  be  distinguished  due  to 
the  boundary  truncation  error  associated  with  the  n=ll  mesh. 


Figures  5-1,  5-2,  and  5-3  show  the  computed  axisymmetric  (p  = e=  1) 
pressure  profiles  for  tt\e  three  material  combinations.  The  limiting 
cases  of  a = 0 and  a = <”  correspond  to  the  Hertz  solutions  for  an 
infinitely  thick  layer  and  a vanishing  layer.  The  variation  of  these 
pressure  profili’s  with  u represents  the  variation  with  layer  thickness 
of  the  load  ri'cpiired  to  maintain  a fixed  contact  zone  radius  a.  When 
each  profile  is  normalized  by  its  center  value,  the  profile  shapes 
shown  in  Figures  5-4,  5-5,  and  5-6  are  obtained.  Those  in  Figure  5-4 
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for  material  combination  (i)  appear  to  approach  a parabolic  shape  with 
increasing  a.  For  material  combination  (iii)  in  Figure  5-6,  very 
little  change  from  the  Hertz  profile  shape  is  evident.  Tliis  is  because’ 
the  relative  surface  displacement  is  dominated  by  the  low  modulus 
homogeneous  solid.  The  profile  shapes  in  Figure  5-5  for  material 
combination  (ii)  show  a substantial  departure  from  the  Hertz  ellipse. 

These  shapes,  which  show  a peak  pressure  near  the  edge  of  contact,  also  , 

appear  in  the  results  of  Chen  and  Engel  [4]  and  of  Gupta  and  Walowit  |61. 

Through  interpolation  of  the  arrays  of  values  of  e,  a /h,  a/a  , b/a  , 

o o o o o , 

d/d^,  and  P / computed  for  the  three  material  combinations  at  the  ' 

above  set  of  values  of  a and  p,  curves  are  generated  to  display  a/a  , 

b/a  , d/d  , and  p /p°  as  functions  of  a /h  at  fixed  values  of  c 

o o o o , 

between  1 and  10.  Figures  5-7,  5-8,  and  5-9  show  the  a/a^  curves  for 

the  three  material  combinations;  Figures  5-10,  5-11,  and  5-12  show  the 

b/a^  curves;  Figures  5-13,  5-14,  and  5-15  show  the  d/d^  curves;  and 

Figures  5-16,  5-17,  and  5-18  show  the  p°/p°  curves.  Since  a , d , and 
o o o o 

p^  defined  by  Equations  3.27,  3.25,  and  3.26  are  constants  that  depend 

on  the  load  W,  these  curves  directly  display  the  dependence  of  a,  b, 

d,  and  p°  on  h and  e. 
o 

The  range  of  a^/h  from  zero  to  the  value  at  which  a given  curve 

terminates  covers  the  range  from  a = 0 to  a = 5 for  that  curve. 

The  value  of  the  starting  point  of  each  curve  (a  /h  = a = 0)  is  obtained 

o 

from  the  Hertz  solution  for  a layer  of  infinite  thickness.  In  the 
limit  a^/h  ->•  “>  which  means  a vanishing  layer,  material  combination  (i) 
is  representing  a contact  between  two  rigid  solids;  for  any  curvature 

ratio  the  contact  zone  is  becoming  a point  with  major  and  minor  half-  , 

widths  and  approach  vanishing  and  the  central  pressure  becoming  infinite. 

In  the  same  limit  for  material  combinations  (ii)  and  (Hi),  the  layered  ^ 

I 

kolld  Is  becoming  a homogeneous  solid  of  the  elastic  substrate  material. 

Th4’  values  of  a/a  , b/a  , d/d  , and  p'^/p^  from  the  Hertz  solution  for  1 

o o o ' o 

' .uhstr.iti*  material  are  sliown  as  limiting  values  approaclied  by  the  j 

..  ^ '.•r  m.iterial  combinations  (ii)  and  (iii).  ' I 

T 
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Elementary  physical  reasoning  explains  most  of  the  trends  shown  by 
these  curves.  For  material  combination  (i)  at  a given  curvature  ratio, 
a decreasing  layer  thickness  reflects  a general  stiffening  of  the 
system  since  the  substrate  is  rigid.  Thus,  at  fixed  load,  tlie 
approach  decreases  (Figure  5-13);  the  major  and  minor  h.ilfwidths 
decrease  (Figures  5-7  and  5-10)  causing  the  contact  area  to  decrease 
and  the  central  pressure  to  increase  (Figure  5-16).  The  opposite 
trend  is  shown  by  material  combinations  (ii)  and  (iii).  The  layer  is 
of  higher  modulus  than  the  substrate  so  the  system  becomes  less 
stiff  with  decreasing  layer  thickness.  This  is  consistent  with  an 
increasing  approach  (Figures  5-14  and  5-15)  , increasing  major  and  minor 
halfwidths  (Figures  5-8,  5-9,  5-11,  and  5-12)  and,  thus,  an  increasing 
contact  area  and  a decreasing  central  pressure  (Figures  5-17  and  5-18) . 
The  approach,  major  and  minor  halfwidths,  and  central  pressure  are 
seen  to  be  less  sensitive  to  layer  thickness  for  material  combination 
(iii)  than  they  are  for  material  combination  (ii).  Since  material 
combination  (iii)  contains  a low  modulus  homogeneous  solid,  a given 
change  in  layer  thickness  makes  less  of  a difference  in  the  overall 
stiffness  of  the  system  than  it  does  for  material  combination  (ii) 
where  the  homogeneous  solid  is  rigid.  In  other  words,  the  relative 
surface  displacement  for  material  combination  (iii)  is  dominated  by 
the  low  modulus  homogeneous  solid,  while  the  relative  surface  displace- 
ment for  material  combination  (ii)  consists  only  of  the  surface  dis- 
placement of  the  layered  solid. 

Figures  5-19,  5-20,  and  5-21  show  curves  of  the  ellipticity  ratio 

(a/b)  versus  a /h  for  the  three  material  combinations.  These  curves 
o 

are  simply  obtained  as  the  ratio  of  values  given  by  the  a/a  and  b/a 

curves.  For  a given  curvature  ratio,  material  combination  (i)  shows  a 

steady  decrease  in  a/b  with  decreasing  layer  thickness;  the  contact 

zone  tends  to  become  more  circular  with  decreasing  layer  thickness. 

This  reflects  the  tendency  of  the  low  modulus  layer  material  to  be 

displaced  in  a direction  parallel  totheminor  axis  of  the  rigid 

elliptical  paraboloid  that  is  penetrating  the  layer.  For  material 

combinations  (il)  and  (iii),  a/b  has  tlie  same  value  for  an  infinitely 

thick  layer  (a  /li  = 0)  as  for  a vanisliing  layer  (a  /h  = “)  since 
o o 
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a/b  depends  only  on  curvature  ratio  in  the  Hertz  solution.  Between 
the  two  extremes,  a/b  increases  to  a maximum  value  and  then  decreases 
with  decreasing  layer  thickness.  Material  combination  (iii)  shows 
very  little  rise  in  a/b  as  expected  from  its  relative  insensitivity 
to  layer  thickness.  The  rise  is  prominent  for  material  combination 
(ii)  with  a/b  at  a particular  e reaching  a maximum  value  that  is  less 
than  the  value  of  e but  at  lower  e very  close  to  it.  As  the  maximum 
value  of  a/b  is  approached,  the  contact  zone  boundary  contour  is 
tending  to  coinc ide  with  a curve  bounding  a cross-section  of  the  rigid 
paraboloid . 


As  a simple  example  to  demonstrate  tlie  application  of  these  results 
to  an  engineering  problem,  consider  the  loading  of  an  elastomer  surface 
layer  (Vj^  = 0.5)  by  the  edge  of  a crowned  disc  as  shown  in  Figure  5-22. 

The  undeformed  separation  profile  near  the  contact  zone  will  be  an 

r 

elliptical  paraboloid  of  principal  radii  of  curvature  R and  R where 

X y 

R is  the  crown  radius  and  R is  the  disc  radius.  Let  the  elastomer  have 
X y 

a low  enough  elastic  modulus  that  the  substrate  it  is  bonded  to  and 
the  disc  can  be  regarded  as  rigid.  Then  the  data  for  material  combina- 
tion (i)  is  applicable.  Consider  the  following  case: 


Disc  Radius: 
Crown  Radius: 
Elastomer : 

Layer  Thickness: 
Load : 


R =2  inches 

y 

R =6  inches 

X 

= 1000  Ib/in.^ 
h = 0.10  inches 
W = 20  lb 


One  wishes  to  determine  the  major  and  minor  halfwidths  a,  b,  approach 
d,  and  central  pressure  p°.  First,  from  Equations  3.27,  3.25,  and 
3.26; 

a = 0.407  inches 
o 

d = 0.028  inches 
o 

p°  = 57.6  Ib/in.^ 


j 

I 

With  e = 3 and  a /h  = 4.07,  a/a  , b/a  , d/d  , and  p /p^  can  be  read 
o o o o o 

j directly  from  Figures  5-7,  5-10,  5-13,  and  5-16,  respectively.  One 

1 finds:  J 

( a/a  = 0.745  * 

o 

b/a  = 0.447 
o 

j d/d  = 0.25 

I o 

p°/p°  = 3.9 

! 

I with  the  result: 

i a = 0.303  Inches  i 

I i 

b = 0. 182  Inches 
, d = 0.007  inches 

I 

I p°  = 220  Ib/in.^ 
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PART  6 

SUMMARY  AND  CONCLUSIONS 

In  this  study,  a numerical  method  is  developed  to  compute  the  pressuri' 
distribution  and  normal  approach  in  a generalized  elliptical  contact 
between  layered  solids.  The  computed  quantities  are  obtained  as  an 
approximate  numerical  solution  to  an  integral  equation  formulated  for 
the  most  general  case  of  Hertz's  assumptions,  that  is,  the  frictionless 
contact  between  arbitrary  surfaces  whose  undeformed  normal  separation 
can  be  approximated  by  the  separation  between  an  elliptical  paraboloid 
and  the  tangent  plane  at  its  vertex.  The  method  is  applied  to  the  con- 
tact between  a homogeneous  solid  and  a layered  solid  consisting  of  an 
isotropic  surface  layer  in  perfect  adhesion  to  an  isotropic  substrate. 

While  plane  strain  and  axisymmetric  solutions  for  this  case  have  been 
reported  in  the  literature  [4-8],  no  solutions  have  been  reported  for 
the  general  three-dimensional  problem  dealt  with  here 

With  the  kernel  function  given  by  a Hankel  transform  inversion  integral, 
the  integral  equation  must  be  treated  by  some  approximate  method  to 
obtain  solutions  in  a useful  numerical  form.  The  method  used  here 
combines  a discretized  representation  of  the  unknown  pressure  distri- 
bution with  an  essentially  exact  numerical  evaluation  of  the  kernel 
function.  Values  of  pressure  are  defined  at  points  on  a rectilinear 
grid  representing  the  contact  zone  and  a linear  function  approximates 
the  continuous  distribution  of  pressure  between  adjacent  points.  By 
expressing  the  integral  equation  at  each  grid  point  for  this  discretized 
distribution,  one  arrivts  at  a system  of  simultaneous  linear  equations 
in  the  grid  point  pressures  with  coefficients  formed  from  integrals  of 
the  kernel  over  right  triangular  regions  between  adjacent  points.  The 
approach  and  the  ratio  of  the  principal  radii  of  curvature  of  the  parabolic 
separation  also  appear  as  unknowns. 

The  key  numerical  task  in  this  analysis  is  the  evaluation  of  a large 
array  of  integrals  of  the  kernel  function  required  to  form  the  coefficients 
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of  the  unknown  pressures.  Integration  in  three  variables,  the  two  space- 
variables  of  the  Integral  equation  and  the  Hankel  transform  variable  i-a, 
is  performed  by  Gaussian  quadrature  with  a resultant  accuracy  of  around 
six  decimal  places.  Once  the  coefficients  are  computed,  the  system 
of  equations  for  the  grid  point  pressures,  approach,  and  curvature 
ratio  is  solved  by  Gaussian  elimination.  In  physical  terms,  the 
solution  permits  one  to  determine  the  major  and  minor  halfwidths, 
pressures,  and  approach  for  a given  load,  layer  thickness,  and  curva- 
ture ratio.  The  dimensionless  load,  obtained  as  the  integral  of  the 
solution  pressure  distribution,  is  equivalent  to  a dimensionless 
major  halfwidth. 

In  setting  up  the  contact  zone  grid  for  a particular  problem,  the 
major  to  minor  halfwidth  ratio  or  ellipticity  ratio  is  defined.  The 
shape  of  the  contact  zone  boundary  contour  between  the  vertices  at 
the  major  and  minor  halfwidths  must  also  be  defined  even  though  the 
true  shape  is  an  unknown  in  the  problem.  In  general,  the  true  shape, 
which  is  an  ellipse  in  the  Hertz  solution  for  homogeneous  solids, 
must  be  determined  iteratively  as  the  shape  along  which  the  solution 
pressure  profile  vanishes.  With  a grid  of  ten  uniform  divisions  along 
the  major  and  minor  halfwidths,  a discrete  representation  of  an 
elliptical  contour  gave  satisfactory  results  for  all  cases  considered  in 
the  study. 

Accuracy  of  the  numerical  method  with  the  above-mentioned  grid  was 
checked  at  limiting  cases  by  comparison  of  predicted  results  with 
available  solutions.  Comparisons  Include  the  Hertz  solution  for  homo- 
geneous solids  [1-3]  and  numerical  results  of  Chen  and  Engel  [4]  for 
axlsymmetric  layered  solids.  The  major  and  minor  halfwidths,  central 
pressure,  and  approach  computed  by  the  present  method  agree  to  within 
one  percent  with  the  comparison  values. 

Numerical  results  generated  over  a range  of  ellipticity  ratios  and 
major  halfwidth  to  layer  thickness  ratios  are  shown  for  three  material 
combinations  as  curves  that  give  the  dimensionless  major  and  minor 
halfwidths,  approach,  and  central  pressure  as  a function  of  layer 
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thickness  at  fixed  curvature  ratio.  The  load  appears  in  scale  factors 
for  the  dimensionless  variables.  At  a given  load  and  curvature  ratio, 
the  ellipticity  ratio  of  the  contact  zone  is  shown  to  decrease  with 
decreasing  layer  thickness  for  a low  modulus  layer  bonded  to  a high 
modulus  substrate  and  increase  with  decreasing  layer  thickness  for  a 
high  modulus  layer  bonded  to  a low  modulus  substrate.  Pressure  profile 
shapes  that  differ  substantially  from  the  ellipsoid  of  the  Hertz  solu- 
tion appear  in  the  case  of  a high  modulus  layer  bonded  to  a low  modulus 
substrate . 

The  same  method  can  be  used  to  compute  the  subsurface  stresses  in  a 
layered  halfspace  under  an  arbitrary  surface  loading.  The  computation 
of  subsurface  stresses  under  a loading  consisting  of  a computed  con- 
tact pressure  distribution  would  be  of  value  as  a means  of  determining 
safe  operating  loads  for  surface-coated  bearing  elements  in  concentrated 
contact.  In  addition,  the  method  can  easily  be  applied  to  other  contact 
configurations  such  as  the  contact  between  two  layered  solids,  layered 
solids  with  interface  conditions  that  involve  slip  and  friction,  and 
solids  that  contain  more  than  one  layer.  The  inversion  integral  for 
the  kernel  function  would  have  a different  form  in  each  case;  everything 
else  would  be  the  same. 
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ANALYTICAL  EVALUATION  OF  BASIC  INTEGRALS 


The  Integral  functions  defined  by  Equation  4.2.17: 
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Ui 

dudv 


(1.2) 


are  evaluated  by  carrying  out  the  double  integration  in  polar  coordinates 
(r,9)  with: 

u = rcosQ 
V = rsin9 


r = 


u +v 


2 


9 = tan  ^(v/u) 

and  the  differential  element  dudv  replaced  with  rdrd9.  Figure  I-l 
shows  the  types  of  regions  which  must  be  considered  with  (a)  applying 
to  Equation  I.l,  (b)  and  (c)  to  Equation  1.2  for  a slope  of  -t , and 
(d)  to  Equation  1.2  for  a slope  of  +o . In  polar  coordinates,  the 
Integrals  over  these  regions  take  the  following  form: 
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s < pt 


Tr/2  s/ (cose+psinS) 
(s,t,-p)  = ; f 

o o 


ir-tan  [t/(pt-s)]  t/sin0  " \ / 

/ f ]rcos9\  drd0 

Tr/2  s/(cos9+psin9)  ) . 

/rsin0\ 


tan  [t/(s+pt)]  s/ (cos0-psin0) 
(s,t,+p)  = / f 

o o 


tan  [t/(s+pt)]  o 


t/sinO  j \ / 

f I /rcos6>  drd0 


rsin9 


Next,  the  straightforward  r integration  is  carried  out  and  the  remaining 
0 integrals  are  adjusted  where  necessary  to  give  all  of  the  lower 
limits  a value  of  zero.  This  results  in: 
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The  above  expressions  contain  the  integrals  of  five  basic  functions  of 
9,  whicri  can  be  evaluated  by  elementary  techniques  to  yield  the 
following : 
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where 


'I  + p' 


(1.5) 


(1.6) 


a = 
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When  the  appropriate  limits  of  integration  in  Equations  1.3  through  1.6 
are  substituted  for  'p  in  these  integrals,  the  final  closed  form 
expressions  for  Equations  I.l  and  1.2  are  obtained.  Defining  the 
quantities: 

R®  = 


97 


I 


APPENDIX  II 

SURFACE  DISPLACEMENT  DUE  TO  AXISYMMETRIC  ELLIPSOIDAL 
PRESSURE  DISTRIBUTION 


Equation  2.12  with  Equation  2.10  glve.s  the  surface  displacement  u^(p) 
of  a layered  halfspace  due  to  an  arbitrary  axlsymmetr Ic  pressure  dis- 
tribution p(p)  where  p » r/h  Is  a dimensionless  radial  coordinate. 

Let  the  pressure  be  distributed  over  a region  of  radius  a and  define 
the  dimensionless  radial  coordinate  > = r/a.  Written  in  terms  of  / , 
Equations  2.12  and  2.10  appear  as: 

2h(l-v 

u (A)  = = / (g  (fj.)-g  (u)]P(w)J  (aA(..)d!.- 

Ef  o 2 4 o 


P(aj)  = / Ap(A)J  (ciAo.)dA 

o 

Now  consider  the  ellipsoidal  pressure  distribution: 


for  which  P('i>)  can  be  shown  to  take  the  form: 

. 2 o sinofij  - 0K.jCosai., 

P(cu)  = T p — r 

permitting  to  be  determined  from: 


u^fA) 


2ap''(l-v^^)  "> 

/ (gT(“)-g/ ('^)  ] J (aAiu)  (sinou.j 

2 4 o 


E 


fiwcosat.)) 


dill 

(aw) 


To  write  this  expression  in  terms  of  the  dimensionless  pressure  de- 
fined in  Equation  1.9,  define  the  dimensionless  displacement: 


which  is  the  scaling  used  implicitly  in  deriving  Equation  3.10.  With 
2R  (1-v  ^p» 

r - 
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APPENDIX  HI 


1 
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DESCRIPTION  OF  COMPUTER  PROGRAM  "HI.IPRG" 


This  program  computes  and  writes  into  a direct  access  file  the  arrays 
of  integral  functions  Hj , 1^,  H*,  I*,  and  J*  according  to  the 

analysis  given  in  Section  4.2.2.  The  program  uses  the  CDC  checkpoint/ 
restart  facility  to  permit  the  £!xecution  of  longer  runs  in  a series 
of  job  steps.  The  parameter  TCHECK,  defined  in  statement  No.  10  of 
HIJPRG,  sets  the  maximum  number  of  central  processor  seconds  for  a single 
job  step.  In  this  version,  TCHECK  = 480.  The  system  clock  reading 
is  compared  with  TCHECK  at  various  points  throughout  HI.TPRG.  If 
TCHECK  is  exceeded,  a checkpoint  dump  is  taken  and  execution  stops. 
Execution  continues  from  that  point  when  the  restart  job  is  submitted. 
The  user  should  consult  the  CDC  FORTRAN  and  NOS  Operating  System  'Manuals 
for  specific  instructions  on  direct  access  file  handling  and  the  check- 
point/restart facility. 


Input  Description 

Punched  card  input  is  submitted  as  follows; 


Card  1 FORMAT  (110) 

ICASE  = Integer  identification  number  for  direct  access 

file  containing  output  data.  For  ICASE  = 0,  output 
is  not  written  on  file. 


Card  2 F0RMAT  (5F10.2) 

P0IS1  = Vj 
P0IS2  = v^ 

E1E2  = Ej/E^ 

Card  3 F0RMAT  (2F10.2,  110) 

AH  = Qi  = a/h 

RH0  = p = a/b 

NDIV  * r - 1 (NDIV  cannot  exceed  10  in  present  versions  of 
HIJPRG  and  PS0LV) 
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Output  Descript  ion 

Printer  output  is  given  for  the  computed  arrays  of  integral  functions: 

XH(I.,I)  = [l^((I-l)Ar.  (.1-1), \nl 

XI(I,.I)  = I|  [(r-DAA,  (,J-l)Ari| 

xj(T..r)  = .I^{(r-])AC.  (.J-I)\nl 

1=1, . . . ,2(NDIV)  + 1 
,1=1 2(NDIV)  + 1 

XHM(I,.J)  = H*  [(I-l)AC,  (J-l)An,  - O] 

XIM(I,.I)  = I*  [(I-l)AC.  (J-l)An.  - ol 

X.IM(I,,J)  = J*  [(I-l)A^,  (J-l)An.  - o] 

i=l,...,4(NDIV)  + 1 
.J=l 2(NDIV)  + 1 

XHP(I,,I)  = H*  ((I-l)Af,,  (J-l)An,  + p1 

XIP(I,.I)  = I*  f(I-l)Ac,  (J-l)An,  + p] 

X.IP(I.I)  = J*  {(I-l)A^,  (,I-l)An,  + p] 

1=1, .. . ,2(NDIV)  + 1 
.1=1 VDIV  + 1 

In  the  output  listing  for  each  array,  "I"  designates  the  number  of 
a line  or  group  of  lines.  The  quantities  appear  sequentially  in  "I" 
in  each  line  or  group  of  lines. 
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DESCRIPTION  OF  COMPUTER  PROGRAM  "PSOLV" 


This  program  reads  the  direct  access  file  created  by  HIJPRG  that 
contains  the  arrays  of  I^,  H*.  I*,  J*  and  completes  the  computa- 
tions described  in  Part  4.  First,  with  read-in  values  of  and 

it  computes  the  arrays  of  H,  I,  .1,  H*,  I*,  and  J*  according  to 
Equations  4.2.15  and  4.2.19  from  the  arrays  of  H^,  1^^,  H*,  I*,  and 

J*  given  by  the  file.  Then,  with  a read-in  contact  zone  boundary 
description,  the  computations  contained  in  Equations  4.2.13,  4.2.6, 
4.2.7,  4.2.5,  and  4.7  and  Table  4-1  are  performed  to  generate  the 
matrix  of  influence  coefficients  C,.  , . , . In  the  notation  of  Equation 
4.14,  the  matrix  inversion  subroutine  MATIN  is  used  to  compute  the 
contact  pressure  distribution  $ .,  approach  P , curvature  ratio  e, 

and,  through  Equation  4.17,  load  W.  The  quantities  a /h,  a/a  , b/a  , 

o o o 

d/d^,  and  p/p  are  computed  from  Equations  3.28  through  3.32. 


The  user  should  consult  the  CDC  NOS  Operating  System  Manuals  for 
specific  instructions  on  handling  the  direct  access  file  created  by 
HIJPRG. 


Input  Description 

Punched  card  input  is  submitted  as  follows: 

Caxd_l  F0RMAT  (110,  2F10.3) 

JCASE  = Integer  identification  number  for  direct  access 

file  from  HIJPRG.  This  must  be  the  same  as  the 

value  of  ICASE  specified  in  executing  HIJPRG. 

E1E3  = E^/E^ 

P^IS3  = 

Card  2 FORMAT  (110) 


ITER 


Integer  identification  number  used  to  identify 
particular  selection  of  contact  zone  boundary. 
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Card  3 F0 RMAT  (1615) 

(MBDY(I),  1=1,  MDIV) 

NBDY(I)  = Number  of  grid  point  along  minor  halfwidth  at 
which  contact  zone  terminates  for  Ith  division 
along  major  halfwidth. 

rr(l  ) = NBOYC  1 ) 

m(I)  = NBOYfr-l)  1 = 2 NDIV  + ! 

Restrictions; 

NBDY(l)  = NDIV  + 1 

NBDY(I)  ^ NBDY(I-l)  , 1=2,..., NDIV 


Card  4 FORMAT  (110) 

NPZERO  = Total  number  of  grid  points  along  boundary  at 
which  zero  pressure  is  to  be  prescribed  (not 
including  points  at  ends  of  major  and  minor  half 
widths  where  pressure  is  automatically  set  to 
zero) . 


Card  _5 

(MZF.ROd)  , 
MZERO(I)  = 

Card  6 

(NZERO(I) , 
NZERO(I)  = 


F0R.MAT  (2014) 

1=1,  NPZERO) 

Vertical  grid  line  (position  along  major  half- 
width) locating  Ith  zero  pressure  boundary  point 

FORMAT  (2014) 

1=1,  NPZERO) 

Horizontal  grid  line  (position  along  minor  half- 
width) locating  Ith  zero  pressure  boundary  point 


Output  Description 

Printer  output  is  given  for  the  following  quantities: 

XH(r,J)  = Hf(I-l)Af,,  (J-l)Anl 

XI(I,J)  - I[(I-l)Af;.  (J-l)Anl 

X.J(I,J)  « J[(I-1)AC,  (J-l)Anl 

I-l 2(NDIV)  + 1 

J-1 2(NDIV)  + 1 


. f 
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XHM(I..J)  = H*| 

(i-DAi;. 

(.1-1  )Ar), 

- p 

XIM(I,.I)  - I*((I-l)Af,, 

(.J-l)An. 

- p 

X.IM(I,.J)  = ,1*1 

(I-l)Ar, 

(.J-l)An, 

- 0 

1 = 1 4(NDIV)  + 1 

.1=1  ....  ,2(NDIV)  + 1 


XHP(I,  J)  = H*[a-11AC,  (J-l)Ar,,  + (j  ] 

XIP(I,J)  = I*[(I-1)AC,  (J-l)An,  + d] 

XJP(I,,I)  = .J*[(I-1)AC,  (J-llAri,  + o\ 

1 = 1 ... . ,2(NDIV)  + 1 
J=l, . . . ,NDIV  + 1 

In  the  output  listing  for  each  array,  "I"  designates  the  number  of  a 
line  or  group  of  lines.  The  quantities  appear  sequentially  in  "J"  in 
each  line  or  group  of  lines. 

Computed  Contact  Parameters: 

kXRY  = £ 

GAM  = r 

PHIO  = (t)° 

WBAR  = W 

ACHH  = a /h 
o 

AACH  = a /a 

o 

BACH  = b/a 

o 

DDCH  = d/d 

o 

POPOCH  = p°/p° 

O 

Hertz  solution  for  infinite  layer  thickness  with  same  values  of  p,  Ej/E^, 


and  v^: 
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RXRYH  = c 
GAMH  = r 
PHIOH  = 4,° 

WBARH  = W 

M=l, . . . ,n 
N=1 m(M) 

Pressure  Distribution  for  Hertz  Solution: 

PHIH(XI(M),ETA(N))  = * „ M=l,...,n 

M,N 

N=l,...,m(M) 

The  pressure  distributions  are  listed  as  two-dimensional  arrays.  The 
line  number  is  designated  by  M and  the  values  on  each  line  are  given 
sequentially  in  N. 


Pressure  Distribution: 

PHI(XI(M),ETA(N))  = 4^ 


APPENDIX  V 


LISTING  OF  COMPUTER  PROGRAM  "PSOLV" 
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